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ABSTRACT 


This  report  deals  with  the  specification  of  the  angular  position 
of  objects  by  measurement  of  the  angles  of  arrival  of  light  beams;  and 
with  the  measurement  uncertainties  due  to  atmospheric  turbulence.  A 
survey  of  the  literature  is  presented  showing  two  approaches  to  angle 
of  arrival,  one  for  large  aperture  receivers  and  one  for  a  single 
small  aperture  or  a  pair  of  pinholes.  Pertinent  defining  equations 
are  presented  for  both  cases  and  values  are  calculated  for  arrival 
angle  mean  square  measured  with  a  large  aperture  and  correlation  of 
arrival  angles  for  small  apertures.  A  synonymity  between  the  small 
aperture  correlation  function  and  large  aperture  mean  square  angle  of 
arrival  is  presented  and  used  to  present  further  large  aperture  mean 
square  values.  This  provides  a  significant  simplification  in  calcu¬ 
lation  and  measurement  capabilities.  Finally  the  regions  of  accuracy 
are  presented  for  the  Rytov  approximation  which  forms  the  basis  for 
the  numeric  calculations. 
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INTRODUCTION 


When  optical  methods  are  used  to  observe  an  object  through  the 
atmosphere,  the  turbulence  in  the  atmosphere  may  cause  the  apparent 
position  of  the  object  to  fluctuate.  This  motion  is  caused  by  turbu¬ 
lence-induced  fluctuations  in  the  electromagnetic  phase  of  the  received 
wave.  The  received  rays  appear  to  change  direction,  so  that  the  term 
angle  of  arrival  fluctuations  is  used  to  describe  this  phenomenon. 
Knowledge  of  the  angle  of  arrival  fluctuations  or  its  statistical 
properties  is  obviously  necessary  for  analysis  of  systems  specifying 
the  accurate  location  of  terrestrial  or  extraterrestrial  object  by 
optical  means.  It  is  also  conceivable  that,  since  angle  of  arrival 
fluctuations  are  induced  by  atmospheric  turbulence,  arrival  angle 
statistics  could  be  used  as  a  convenient  tool  for  study  of  turbulence 
characteristics. 

The  literature  on  angle  of  arrival  is  scattered  and  offers  a  wide 
range  of  views,  approaches,  and  restrictions.  For  example  there  are 
two  general  approaches  to  arrival  angle.  The  first  considers  the  light 
detected  by  one  or  several  minute  apertures  where  phase  is  linear  and 
amplitude  is  constant  across  the  aperture.  The  quantities  of  interest 
for  this  approach  are  the  arrival  angle  variance  and  correlation 
function,  the  mean  being  easily  demonstrated  to  be  zero.  Alternatively 
there  is  consideration  given  to  a  single  large  aperture  where  the  phase 
front  is  crinkled  and  the  amplitude  is  not  at  all  constant  across  the 
aperture.  Generally  the  literature  considers  only  the  arrival  angle 
variance  for  such  a  situation.  For  either  approach  the  derivations 
based  on  the  Rytov,  the  geometrical  optics  or  other  approximation  may 
be  gi\en. 

It  is  desirable  to  be  able  to  predict  the  angle  of  arrival 
observed  by  both  large  and  small  aperture  methods  in  terms  of  the  propa¬ 
gation  range,  receiver  and  aperture  characteristics,  wavelength,  and 
turbulence  parameters;  outer  and  inner  scale  and  structure  parameter. 

Some  cases  have  been  considered,  but  most  are  valid  only  for  restricted 
ranges  of  the  variables  described  earlier  and  all  use  several  approxi¬ 
mations. 

The  original  object  of  this  report  was  to  examine  the  literature, 
to  categorize  these  different  situations  and  approximations,  and  to 
extend  the  results  where  feasible.  This  has  been  done.  Basic  extensions 
involve  the  computation  of  arrival  angle  functions  using  an  index 
spectrum  which  accounts  for  the  turbulence  inner  and  outer  scales. 

In  the  process  of  performing  these  calculations,  a  synonymity 
between  a  small  aperture  correlation  function  and  large  aperture  variance 
function  was  discovered.  This  synonymity  provides  great  simplification 
in  the  calculation  and  measurement  of  these  two  quantities.  The 
synonymity  is  presented  and  derived. 
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In  the  balance  of  the  report  there  is  first  a  general  discussion 
of  angle  of  arrival  in  Chapter  II.  This  is  followed  by  analytic 
expressions  for  angle  of  arrival  statistical  functions  for  first  the 
small  aperture  approach  (Chapter  III)  and  then  the  large  aperture 
approach  (Chapter  IV).  Numerical  results  are  presented  in  Chapter  V 
in  graphical  form  for  a  wide  variety  of  cases  along  with  a  discussion 
of  the  limits  of  appl i cabi 1 i ty  of  the  Rytov  approximation.  Chapter  VI 
deals  with  the  synonymity  of  the  large  and  small  aperture  approaches. 
Finally,  the  summary  appears  in  Chapter  VII. 


II.  GENERAL  DISCUSSION  OF  ANGLE  OF  ARRIVAL 

In  this  section  we  introduce  in  detail  the  various  concepts  and 
quantities  associated  with  the  two  major  approaches  to  angle  of 
arrival.  First  the  general  measurement  situation  is  considered.  The 
small  and  large  aperture  situations  are  discussed  to  the  extent  that 
typical  formulae  are  presented. 

As  generally  used,  the  term  angle  of  arrival  is  concerned  with 
light  emitted  from  an  object  and  detected  by  a  receiver.  Pragmatically 
what  is  desired  is  a  listing  of  the  angular  coordinates  of  the  object 
with  respect  to  some  fixed  reference  system.  In  simple  minded  terms 
this  can  be  obtained  by  imaging  the  object  either  with  a  lens  (finite 
aperture)  or  a  pinhole  camera  (very  small  aperture)  and  taking  a  line 
from  the  object  through  the  center  of  the  input  aperture.  The  angular 
coordinates  of  the  line  then  define  the  direction  of  the  object.  If 
there  are  no  intervening  refractive  index  fluctuations  then  the  image 
of  a  point  object  may  be  a  diffraction  limited  spot  with  a  well  defined 
center.  The  angular  coordinates  determined  by  the  spot  center  can  be 
very  precisely  chcsen  in  such  a  case. 

Complications  arise  because  of  random  index  fluctuations  which 
cause  the  image  to  move  around  and/or  to  blur,  depending  on  the  physical 
details  of  the  receiver.  Random  image  motion  is  then  interpreted  as 
random  fluctuations  in  the  observed  angle  of  arrival.  Further  fluctu¬ 
ations  in  the  shape  of  the  received  image  of  a  point  object  also  cause 
fluctuations  in  the  "center  of  gravity"  of  the  image  and  may  be  in¬ 
terpreted  as  fluctuations  of  the  angle  of  arrival.  It  is  the  statistical 
description  ov  these  random  fluctuations  in  arrival  angle  caused  by 
refractive  index  fluctuations  that  form  the  basic  subject  of  this  report. 

It  might  be  noted  that  the  arrival  angles  measured  with  an  optical 
instrument  might,  even  in  the  absence  of  recractive  index  fluctuations, 
not  be  the  true  angular  coordinates  of  the  object.  For  example,  if 
there  were  a  refractive  index  gradient  perpendicular  to  the  light  path, 
the  light  rays  would  be  bent  causing  an  anomolous  indication  of  the 
true  angles.  This  effect  causes  the  phenomenon  of  "looming"  and 
causes  the  sun  to  be  visible  after  it  has  in  actuality  dropped  below 
the  horizon.  Such  effects  will  not  be  considered  in  the  present  report. 
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There  are  various  approaches  to  the  measurement  and  prediction  of 
angle  of  arrival  fluctuations.  The  measurements  are  generally  cate¬ 
gorized  by  the  size  of  the  aperture  used  for  the  receiver.  Perhaps  the 
cleanest  approach  involves  the  use  of  a  very  small  pinhole  sized  aperture, 
as  shown  in  Fig.  1.  The  pinhole  is  chosen  sufficiently  small  so  as  to 


y  DIFFRACTION 


Fig.  1.  Definition  of  small  aperture  angle  of  arrival. 


be  smaller  than  the  spatial  extent  of  any  amplitude  variation,  and  so  as 
to  contain  only  a  flat  wavefront.  The  displacement,  d,  of  the  center  of 
the  diffraction  pattern  in  a  given  direction,  then,  is  related  to  the 
angle  of  arrival,  a,  by  the  expression 


(1)  d  =  Ra 


where  R  is  the  distance  from  pinhole  to  diffraction  pattern.  In  this 
situation,  motion  is  caused  only  by  phase  fluctuations  and  not  by 
amplitude  fluctuations  of  the  incoming  light. 
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As  shown  in  Fig.  1,  two  angles  of  arrival  are  defined;  the  ele¬ 
vation  angle  a,  and  the  azimuth  angle,  &.  One  can  consider  various 
statistical  functions  of  these  angles.  The  simplest  are  the  variances,* 
such  as  <(a-<a>)2>.  (The  angular  brackets  indicate  ensemble  average.) 
There  are  also  the  two  point  arrival  angle  correlation  functions  for 
angles  of  arrival  defined  at  two  separate  pinholes  as  shown  in  Fig.  2. 
Two  such  correlation  functions  exist,  differentiated  by  the  relative 
positions  of  the  pinholes  with  respect  to  the  direction  of  motion  of 
the  spots,  as  shown  in  Fig.  2.  These  will  be  discussed  in  more  detail 
in  the  next  section. 

To  calculate  angle  of  arrival  for  a  pinhole  aperture  the  wave- 
front  tilt  is  determined  using  the  optical  distance,  L.  or  the 
equivalent  electromagnetic  phase  shift,  S  =  k  Lt  along  two  rays  through 
the  turbulent  atmosphere  from  the  source  to  opposite  sides  of  the 
receiver  aperture.  As  shown  in  Fig.  3,  we  have  for 


(2) 


a  ,  Um  j-.toJ_-.kCjd.  =  1L  =lii 
“  n  3y  k  3y 

ry>0 


The  type  of  wave,  i.e.,  plane,  spherical,  etc.,  will  determine  the 
functional  form  of  L.  The  arrival  angle  variance  «* ^  is  then  expressed 
as  the  ensemble  average  of  the  square  of  the  average  product  of  two 
optical  phases. 


(2a) 


^  =  ?  ^2  <S{y])  S(y2)> 


yry2 


(2b) 


(2c) 


^2  3yi93y2  Bs(y2'yl) 


yl=y2 


-1 


yry2 


*The  terms  angle  of  arrival  variance  and  mean  square  value  will  be  used 
interchangeably.  They  are  identical  since  the  mean  is  chosen  to  be 
zero.  Similarily,  the  terms  covariance  and  correlation  function  will 
be  used  interchangeably. 
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DIFFRACTION  PATTERN 


PINHOLES 

(b) 

Fig.  2.  Definition  of  small  aperture  angle  of  arrival 
correlation  functions. 
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INCOMING  OEGRAOEO 

WAVEFRONTS  WAVEFRONT 


Fig.  3.  Example  of  a  small  aperture  angle  of 
arrival  calculation. 


Equation  (2a)  uses  the  standard  form  for  the  phase  correlation  function 


B^y^)  =  <S(y])  S(y2)> 
and  the  associated  structure  function, 


Ds(y2-y,)  -  2  [Bs(0)  -  8£ (y g-y, )] 

written  for  a  homogeneous  atmosphere.  Equations  similar  to  Eq.  (2c) 
will  be  used  subsequently  in  the  report. 

Various  techniques  appear  in  the  literature  for  the  extensions 
and  evaluation  of  Eq.  (2c).  They  use  expressions  for  the  phase 
correlation  or  structure  functions  based  on  either  the  eikonal  equation 
approachp  ,2,3]  of  ray  optics,  or  based  on  the  Rytov  approximation^, 5, 6, 7] 
(which  reduces  to  the  eikonal  equation  approach  for  short  ranges).  They 
also  use  a  variety  of  correlation  functions,  structure  functions  or 
spatial  spectra  to  represent  the  atmosphere. 
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An  apparatus  using  pinhole  apertures  has  one  possible  drawback  as 
a  practical  device  for  angular  coordinate  measurement.  That  is  the 
high  source  signal  level  required  because  of  the  small  receiver  aperture 
size.  However  It  does  provide  a  measure  of  arrival  angle  variance  and 
correlation  functions  and,  therefore,  could  be  useful  under  controlled 
conditions  where  the  source  is  available.  Further  the  associated  ana¬ 
lytical  calculations  of  the  statistical  quantities  are  simplified  by 
being  amplitude  independent. 

The  other  approach  to  arrival  angle  measurement  employs  a  large 
aperture  receiver  with  a  lens  to  image  the  object  as  shown  in  Fig.  4. 


Fig.  4.  Definition  of  large  aperture  angle  of  arrival. 


This  situation  is  more  practical  because  of  its  light-gathering  ability 
and  subsequent  increased  sensitivity.  However,  it  is  more  difficult 
to  work  with  conceptually  because  of  random  spatial  wavefront  and 
amplitude  fluctuations  superimposed  upon  the  no-turbulence  fields. 

Thus,  while  the  wavefront  normal  provided  a  reasonable  indication  of 
arrival  angle  for  a  pinhole  aperture,  there  is  very  apt  to  be  no  unique 
wavefront  normal  for  a  large  aperture.  One  can  instead  use  the  focused 
image,  but  amplitude  and  phase  fluctuations  about  the  mean  wavefront 
of  the  'icoming  light  cause  motion  of  the  image  and  deterioration  of 
its  shape,  thereby  increasing  measurement  uncertainty. 
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The  quantity  generally  calculated  for  the  large  aperture  approach 
is  the  arrival  angle  variance.  To  do  this  the  instantaneous  arrival 
angle  is  defined  using  some  sort  of  optimization  technique  involving 
either  fitting  a  smooth  wavefront  to  the  input  fields[8,9]  or  defining 
the  arrival  angle  to  be  that  which  has  the  maximum  power. [10]  These 
calculations  are  generally  more  complex  than  pinhole  aperture  calculations 
because  of  amplitude  fluctuations  over  the  input  aperture,  and  because  of 
the  multiple  integrations  required. 

As  an  indication  of  the  approach  to  arrival  angle  calculations 
for  large  aperture  situations,  consider  the  following  approximate 
calculation  from  Hufnagel.[8]  Assume  that  to  a  good  degree  of 
approximation  the  wavefront  at  any  instant  can  be  described  by  a  plane 
which  matches  the  actual  wavefront  at  four  points,  the  intersection  of 
the  transverse  coordinate  axis  with  the  aperture  edge.  Then  the  phase 
deviation  along  the  x  axis  will  be 

0)  -  4>{-^,  0'/, 


the  linear  variation  of  the  wavefront  from  the  x  axis  will  be 

[+(§,  0)  -  *(-§,  0) ]/k 


and  the  angle  of  tilt  in  the  x  direction  will  be 


(3a) 

cx  =  [»(£.  0)  -  ♦(-§,  0)]/kD. 

Assuming 

that  a  has  a  zero  mean,  then  its  variance 

X 

(3b) 

<ax2>  =  <U(f,  0)  -  *(-£,  0)>2/k2  D2 

(3c) 

=  Ds(D)/k2  D2 

(3d) 

=  6.88(D/rQ)5/3/k2  D2 

(3e) 

=  6.88/k2  rQ5/3  D1/3 

where  the  phase  structure  function,  Ds(p),  has  been  approximated  by 
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Os(p)  =  6.88(p/r0)5/3 

where  rQ  is  Fried's  transverse  coherence  length  given  by 
rQ  =  (6.88/2.91  k2L  C2)3/5. 

There  are  also  other  approaches  to  angle  of  arrival.  One  is 
an  examination  of  arrival  angle  of  light  from  a  line  source, [11 ,12] 
representing  an  extension  of  the  small  aperture  approach.  Another  is 
a  calculation  of  arrival  angle  variance  using  the  Fokker-Plank 
equation[13]  which  represents  a  unique  approach.  These  techniques 
will  not  be  mentioned  further.  A  few  experimental  measurements  have 
been  described  by  Coulman. [14,15] 

To  summarize,  this  section  has  presented  a  general  discussion  of 
various  aspects  of  angle  of  arrival  work,  indicating  the  sources  of 
uncertainty  considered,  and  discussion  of  the  two  models  for  approaching 
the  subject. 


III.  SMALL  APERTURE  ANGLE  OF  ARRIVAL  CALCULATIONS 

In  this  section  the  exact  expression  for  the  small  aperture 
angle  of  arrival  statistical  quantities  are  derived  and  examined. 

They  include  two  types  of  correlation  functions  and  the  mean  square 
angle  of  arrival.  The  formalism  to  be  presented  follows  closely  that 
presented  by  Tatarski  [5,  section  40]*.  As  mentioned  earlier,  the 
small  aperture  technique  assumes  that  the  aperture  is  small  enough 
that  the  aperture  size  has  no  effect  on  the  angle  of  arrival  indication. 
Thus,  the  phase  difference  across  the  aperture  is  assumed  linear  and 
the  instantaneous  amplitude  is  assumed  to  be  independent  of  position 
in  the  aperture.  Under  these  conditions  the  wavefront  normal  will 
uniquely  define  the  angle  of  arrival  at  any  instant  of  time. 

Since  the  aperture  size  will  not  affect  the  final  result,  the 
optical  path  lengths  to  pairs  of  points  are  considered  and  the  aperture 
sizes  shrink  to  zero.  The  situation  is  shown  in  Fig.  5a.  Using  the 
definition  in  Eq.  (2)  we  have 


*For  equations  taken  from  [5],  the  section  number  and  equation  number 
will  be  given  in  parentheses. 
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Fig.  5.  Geometry  used  to  calculate  small  aperture 
arrival  angle  correlation  functions. 
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(4a) 


where 

S-|  =  phase  at  the  point  (z,0,0) 

52  =  phase  at  the  point  (z.O.r^) 

53  =  phase  at  the  point  (z,0 .y-|-n2/2) 

54  =  phase  at  the  point  (z.O.y^n^). 

These  points  are  shown  in  Fig.  5a. 

The  elevation  angle  of  arrival  correlation  is  defined  as 

(5a)  Ba(yi )  =  <a(z,0,0)  a(z,0»y^)>. 

To  simplify  we  substitute  Eq.  (4a)  and  Eq.  (4b)  into  Eq.  (5a) 
giving 
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(5b) 


B  (y) 

ct  * 


=  limit 
n-|->0 


<(S1 


S2) (S3 


SA)> 


K  T)-j 


where 


nl  =  n2‘ 


Uiing  the  identity  that 

(5c)  (S]  -  S2)(S3  -  S4)  =  1/2[(S1  -  S4)2  +  (S2  -  S3)2  -  (S1  -  S3)2 


-(s2  -  s4)2]. 

Equation  (5b)  can  be  written  in  terms  of  the  phase  structure  function, 

Ds(p). 

(5d)  Bjy^  =  limit  — 2 — 2  ^°s^yl+nl ^  +  °s^Tnl)  ‘  2Ds^yl^‘ 

^-*0  2k  nl 

Expanding  t>s  (y^  ±  n-i)  in  a  Taylor  series  we  have 

(5e)  D$(yi  ±  n])  =  Ds(yi)  ±  ni  D^(y^)  +  1/2  ^  DJ^)  +  ••• 


which,  when  substituted  into  Eq.  (5d) ,  yields  [5,  (40.35)]  written 
as  Eq.  (5f) 


(5f>  W  =  ^r°s  on) 


where  D^(y)  is  the  phase  structure  function  and  the  primes  denote 
derivatives  with  respect  to  the  argument. 

The  azimuth  angle  of  arrival,  6,  shown  in  Fig.  (5b)  is  calculated 
in  a  similar  manner.  Using  the  points  shown  in  Fig.  (5b)  we  have 


(6a) 


8-j  =  limit 

Cl+O 


(S,  -  s2) 
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(6b) 


eg  85  limit 
?2^0 


where 

S-|  ■  phase  at  the  point  (z,  -  0) 

52  s  phase  at  the  point  (z,  +  0) 

53  =  phase  at  the  point  (z,  -  c2/2»  y-|) 

54  =  phase  at  the  point  (z,  +  ?2/2.  y-|) 

and  with  ^ 

(6c)  B^)  =  lim1t  '  2~V~  <^S2  '  S1^S4  "  s3^- 

5,-fl  k  51 

Using  the  identity  in  Eq.  (5c)  we  again  have  the  correlation  in  terms 
of  the  phase  structure  function.  The  only  difference  is  that 


Then  Eq.  (6c)  becomes 

(6e)  Bg(y1 )  =  limit  ^  [Ds  (\l yf  +  cf )  -  D$(y1)]. 

C,^0  k  51 

Since  the  limit  will  force  5  to  be  arbitrarily  small  the  square  root 
in  Eq.  (6e)  can  be  expanded  in  a  power  series. 

(6f)  (J  yl  +  ?1  )  =  yl  +  2y^+  **•  * 
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Using  the  Taylor's  series  expansion  in  Eq.  (5e)  we  have 


(6g) 


(GR)- 


W  Ds(yl}  + 


When  substituted  into  Eq.  (6e)  the  resulting  expression  is  [5,(40.38)], 
written  as  Eq.  (6h) 


(6h) 


Vyi}  =  1,2 


W 


2k 


Comparing  Eqs.  ( 5f)  and  (6h)  we  see  that  the  interpair  displace 
ment  with  respect  to  the  pair  separation  determines  the  final  form  for 
Ba(yi)  and  Bg(y|).  To  calculate  Ba(y|)  the  displacement  was  parallel 
to  the  separation.  On  the  other  hand,  the  displacement  was  transverse 
to  the  separation  when  Bg(yj)  was  calculated. 

The  dependence  on  the  interpair  displacement  when  the  interpair 
direction  is  neither  parallel  or  perpendicular  to  the  interpoint 
direction  has  been  found  explicitly  by  Strobehn  and  Cl  if ford. [7]  For 
sample  points  oriented  as  in  Fig.  6  we  have 


Fig.  6. 


Generalized  geometry  to  calculate  small  aperture 
arrival  angle  correlation  functions. 


(7) 


B£(p) 


~  T  D"(p)  COS^y  +  — T 

2k^  s  2k 


D^(p) 

P 


sin 


Y 


where 

e  =  angle  of  arrival  of  interest 

Y  =  angle  between  interpair  displacement 
and  pair  separation. 


For  y  =  0  and  ir/2  we  have  the  expressions  for  Ba(p)  and  B0(p) 
respectively. 

The  smallness  criterion  for  the  aperture  size  can  now  be 
considered.  Recall  that  it  was  assumed  that  the  aperture  size  had 
no  effect  on  the  angle  of  arrival.  From  Eqs.  (5f)  and  (6h)  we  see  that 
for  Ds(p)a  pz,  the  angle  of  arrival  correlation  'is  independent  of  p. 

It  is  well  known  that  for  the  Kolmogorov  refractive  index  fluctuation 
spectrum,  Ds(p)  has  this  form  for  p  less  than  the  turbulence  inner 
scale.  Thus,  these  results  apply  to  apertures  smaller  than  the  inner 
scale.  For  apertures  larger  than  this,  the  large  aperture  method 
should  be  used. 

In  this  section  the  basic  equations  for  small  aperture  angle  of 
arrival  correlation  functions  have  been  derived.  These  expressions 
will  be  used  in  sections  V  and  VI. 


IV.  LARGE  APERTURE  ANGLE  OF  ARRIVAL  CALCULATIONS 

In  this  section  pertinent  formulas  for  the  large  aperture  angle 
of  arrival  mean  square  value  are  derived  and  discussed.  The  linear 
phase  and  constant  amplitude  assumptions  made  for  the  small  aperture 
angle  of  arrival  calculations  begin  to  break  down  as  the  aperture 
becomes  larger  than  the  inner  scale.  The  phase  front  becomes  more 
distorted  than  can  be  accurately  described  with  a  linear  wave  tilt. 
Secondly,  amplitude  distortions  become  significant.  These  problems 
have  been  considered  in  part  by  various  workers. 

Two  major  approaches  are  used,  both  involving  an  extremum 
technique.  The  first  approach  to  define  the  Instantaneous  arrival 
angle  as  the  direction  for  which  the  instantaneous  power  is  the 
maximum.  The  second  is  to  fit  expression  for  the  instantaneous  wave- 
front  to  a  series  of  polynomial  surfaces  using  a  least  squares  fit. 
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One  of  the  polynomials  is  a  flat  plane,  so  that  the  normal  to  the 
plane  then  is  taken  as  the  direction  of  the  wave.  Both  approaches 
give  the  same  basic  expressions  for  the  instantaneous  angle  of  arrival 
direction.  The  derivations  will  now  be  outlined. 

The  first  approach  to  be  considered  is  that  which  uses  the 
direction  of  maximum  instantaneous  power  to  define  angle  of  arrival [10]. 
The  technique  is  to  split  the  mathematical  expression  for  the  electro¬ 
magnetic  phase  into  two  parts,  one  representing  a  plane  wave  and  the 
other  representing  the  deviations  from  the  plane  wave.  The  direction 
of  the  plane  wave  normal  is  chosen  to  maximize  the  expression  for 
the  power.  It  is  assumed  that  the  deviations  of  the  electromagnetic 
phase  from  the  plane  wave  value  are  small  and  that  amplitude  fluctu¬ 
ations  over  the  input  aperture  are  neglegible.  Thus,  let  the  field 
arriving  at  the  receiving  aperture  be 

EeH.(F)_ 


where  r  =  (x,y)  denotes  the  transverse  position  of  the  input  aperture. 
The  term  <f>(r)  is  equivalent  to  S(x,y,z)  in  the  previous  section.  Then 
the  component  of  the  field  in  the  direction  with  transverse  spatial 
frequency  <  is 


(8)  e(iT)  =  /  dfEei(*(r)  "  K  *  ^ 


and  the  power  in  the  same  direction  is 


(9) 


P(F)  =  i 


I  dFEei(*(r)  *  K  •  r) 


Following  Heidbredder,[10]  we  define  A(r)  the  phase  difference  from 
a  flat  plane. 


A(r)  =  <j>(r)  -  <  •  r. 


Our  first  object  is  to  find  an  expression  for  the  value  of  k 
for  which  P(ic)  will  be  maximum.  We  assume  that  the  phase  difference 
A(r)  from  that  optimum  value  will  be  small.  Then,  rewriting  Eq.  (9), 
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(10a) 

.2 

PGn-H/dFEel4m 

Lo  |  j 

(10b) 

P(ID  “J-  jV  dr  E  cosA(r)J-  +  ^-  -j  /  dr  E  sinA(r)j' 

(10c) 

2  2 

4  h  {/ dF  (i "  ^~) E}  1/ dF  aE  } 

where  the 
respect  to 

r  dependence  of  a  is  implied.  Differentiating  Eq.  (10c)  with 
<x  to  find  the  optimum  value  we  find 

(11a) 

t  ■  0  ■  2  /  dF1  (1  -  r>  E  /  <"-2  It  r- 

A  A 

+  /  drj  aE  /  dr2  f£-  E 

(11b) 


■  2  /  dr,  E  /  dr2  ($2  -  k  •  r2)  x2  E 


-  /  dr  E  /  dr2  Ex2. 


The  aperture  is  assumed  to  be  symmetric  so  that  the  integrals  over  odd 
powers  of  x  and  y  are  zero.  Thus  the  integral  in  the  second  term  and 
the  integral  over  x«y,  in  the  first  term  vanish.  Solving  for  <  then 
gives  x 


(12) 


=  /  dr.Jfr.(r)  x  E 

/  dr  x2  E 
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The  mean  value  of  <x  is  taken  to  be  zero.  This  is  based  on  the  assump¬ 
tion  that  the  atmospheric  turbulence  is  isotropic  so  that  the  mean 
value  of  <}>(r)  is  zero,  i.e.,  fluctuations  of  phase  are  equally  likely. 

The  variance  of  <x  and  the  associated  component  of  the  angle 
of  arrival ,  ax 


(13) 

are  given  by 


(14) 


<(V 


2 

<K  > 
X 


II  drj  dr2  E]  E2  ^(r^  <{>(r2)>  x]x2 
k2|/  dr  x2  E  2 


Identifying  the  phase 


correlation  function 


Bs ( r2“ rl )  =  <i)(ri )  ♦(r2)> 
relating  it  to  the  phase  structure  function 


(15) 


Bs( r2-rl  ^  =  1  C°s(")  "  D,(Vri)] 


s'  2  ’1 


and  again  using  the  fact  that  the  aperture  is  symmetric  gives 


(16) 


2 

"V 


II  dr]  dr2  E1  E2  D^r^r,)  x2 


2k2  /  dr  x2  E 


Equation  (16)  is  generally  evaluated  for  circular  apertures.  To 
obtain  the  appropriate  expression  for  that  case  first  define 


18 


rl =  h 


r2  r2 


x-j  =  r-j  cos  e-| 


x2  =  r2  cos  e2 


so  that  the  denominator  of  Eq.  (16)  is  (taking  the  amplitude  E  to  be 
constant) , 


R  2lT 

ft  _  - 

(17) 

CM 

-iC 

CM 

dr-j  de  r^  cos2e  E 

0  0 

2k2  (,£  )2  E2. 


The  expression  for  the  arrival  angle  variance  is  then 
R  2ir  R  2tt 


M  kll  M  tall  I  —  _  - 

r2  dr-|  de  J  J  r^d^degE-^  cos(e1-e2)D^(>lr^+r2-2r1  r2cos(02-e-| )) 


(18)<aS  =  - 


0  0 


0  0 


2k2  (uR4/4)2  E2 


After  further  defining 
Y  =  ©2_01 

g  =  l/2(e2+©1 ) 


and  performing  the  8  integral,  the  result  is, 

R  R 


(19) 


2  o 

<a  >  =  -  - 

x 


K  K  TT  - - 

|  dr^r2  J  dr2  r2  J  dY  cosy  Ds(Jri'+r2"2rir2  cosy) 


tt  k2(R4/4)2 
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Equation  (19)  will  be  used  in  the  next  section  for  numerical  evaluations. 
Heidbreder  evaluates  Eq.  (10)  with  various  synmetric  aperture  weighting 
functions  W(r])W(r?)  multiplying  the  phase  structure  function  finding 
relatively  small  changes  for  gaussian  weighting  functions  with  half¬ 
widths  the  order  of  the  radius.  He  also  finds  some  small  correlation 
between  the  plane  wave  fluctuations  of  the  phase  front  and  the  residual 
phase  for  a  one  dimensional  aperture  (infinitesimal  in  the  other 
dimension). 

The  second  approach,  due  to  Fried, [9]  in  which  the  expression 
for  the  instantaneous  wavefront  is  fitted  to  a  series  of  polynomials 
will  now  be  considered.  The  polynomials  are  similar  to  the  Zernike 
polynomials,  orthogonal  over  a  circular  area.  The  first  six  are 


they  are  orthonormal  so  that 


(21)  /  dF  K(r>  F„(F)  Fm(f)  -  sm 

where 


(22) 


r  <  R 
r  >  R 
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The  wavefront  4>(r)  is  then  expressed  as  an  infinite  series  of 
these  polynomials: 


*(F)  ■  j,  an  Fn(F)' 


The  second  and  third  polynomials  represent  the  flat  of  the 
wavefront.  Thus  the  spatial  frequency  and  arrival  angle  associated 
with  the  flat  componant  are  given  by 


(23a)  o  =  £  -  grad  {a2  F2(r)  +  a3  F3(F)} 


=  ■1-_  i'a?  x  +  a^  y} 

A  A 

(x  and  y  are  unit  vectors.)  Further,  the  variance  of  the  angle  of 
arrival  in  the  x  direction  is 


(24) 


<a2>  = 
X 


irR4/4 


2 

<a2> 


The  polynomial  series  is  fitted  to  the  wavefront  at  any  instant 
by  choosing  the  coefficients  so  as  to  minimize  the  mean  square 
deviation  between  the  series  and  the  wavefront  proper: 


<25>  I  dr  !♦<»■>  -  1.  an  ^(r)]2  *  0- 


Squaring  the  term  in  the  brackets,  performing  the  indicated  differenti¬ 
ation,  and  solving  for  a2  gives 


(26)  a2  =  /  dr  W(r)  $(r)  Foo(r). 
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Equation  (26)  can  be  shown  to  be  comparable  with  the  expression 
by  Heidbreder.  To  do  that,  use  the  expression  for  F2(7r)  in  Eq.  (20) 
to  give  a£  in  Eq.  (26)  and  substitute  into  Eq.  (24).  After  using 
Eq.  (15)  to  simplify,  the  result  is 


(27) 


2 

<a  >  = 
x 


— U — o  //  dr,  dr2  W(r,)  W(rJ  x,x?  D  (ryrT) 
2k2(,rR4/4)‘-  1  1  1  1  w  s  2  l 


which  coincides  very  nicely  with  Eq.  (18). 

Fried  and  Heidbreder  both  evaluate  the  angle  of  arrival  variance 
assuming  that  the  wave  structure  function  is  approximated  by  the 
phase  stricture  function,  i.e.,  the  range  is  sufficiently  short  that 
amplitude  effects  are  negligible: 


(28)  Dytr)  =  0s(r)  =  6.88(£-)5/3 

where  rQ  is  Fried's  transverse  coherence  length  given  by 


(29)  rQ  *  (6.88/2.91  k2LC^)3/5  . 


The  result  is 


(30a)  <a2>  =  7.064/k2  rQ5/3  D1/3 

(30b)  =  1.026  x  6.88/k2  rQ5/3  D1/3. 


Equation  (30b)  agrees  closely  with  the  approximate  expression 
derived  by  Hufnagel[8]  (Eq.  (3e)). 

To  account  for  amplitude  effects  the  phase  structure  function 
can  be  replaced  with  the  wave  structure  function  and  the  polynomial 
series  now  has  complex  coefficients. [9]  One  drawback  as  explained 
by  Fried  is  that  the  results  are  difficult  to  interpret. 
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Concluding,  we  see  that  the  angle  of  arrival  for  a  large 
aperture  depends  on  the  phase  structure  function  also.  Thus,  these 
results  are  as  general  as  the  approximations  used  to  calculate  the 
phase  structure  function.  A  specific  phase  structure  function  more 
general  than  Eq.  (28)  will  now  be  used  to  calculate  extended  numeric 
resul ts . 


V.  NUMERIC  CALCULATIONS  OF  ANGLE  OF  ARRIVAL 

In  this  section  the  general  expressions  given  previously  for 
arrival  angle  mean  square  and  correlation  functions  are  evaluated. 

The  effects  of  such  parameters  as  outer  and  inner  scale,  range  and 
type  of  wave  are  demonstrated  for  both  large  and  small  aperture  cases. 
The  range  of  parameters  covered  is  such  as  to  give  a  reasonably  complete 
set  of  curves  for  comparison  with  experiments  10. 6y  performed  at  the 
RADC  Laser  Propagation  Range. 

To  proceed  with  the  evaluation,  we  not®  that  the  angle  of  arrival 
statistical  functions  in  Eqs.  (2c),  (5f),  and  (6h)  are  all  in  terms  of 
the  phase  structure  function.  To  evaluate  these  expressions  an  integral 
form  of  Ds(p)  is  used, [5] 


(31)  Ds(p)  *  4„2k2L  |  [1  -  J0(kP)] 

0 

•  $N(<)  <  die 


1  + 


k 

A 


sin 


A 

T 


where 


p  -  separation 
L  =  propagation  range 

$n(k)  =  three  dimensional  refractive  index  fluctuation 
spectrum. 


For  a  spherical  wave,  a  similar  expression  was  derived  by  Carlson 
and  Ishimaru.[16] 
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(32) 


2 


L  r 


Ds(p)  =  sA2  |  [1  -  JQ (,cp ) ]  j 


'klknl^ 


5-  “s  > 


$N  (“4  dn  <cIk, 


Both  of  these  expressions  are  the  result  of  using  the  Rytov  approxi¬ 
mation  and  are  valid  where  the  Rytov  approximation  is  valid.  This 
validity  will  be  discussed  later. 

Both  expressions  (31)  and  (32),  use  the  refractive  index 
fluctuation  spectrum  $«(<).  Several  models  for  $n  have  been  proposed 
and  these  are  discussed  in  Appendix  A.  For  the  calculations  in  this 
report,  the  model  used  is 


(33) 


2  e  m 

$n(k)  =  .033  CN  —2—  2vl1/6 

(<  +  ) 


where  Km  =  5.92/t0 

=  inner  scale  of  the  turbulence 
o 

*0  ■ 

L  =  outer  scale  of  the  turbulence, 
o 


Using  Eq.  (31)  in  the  expressions  derived  earlier  in  Eqs.  (2c), 
(5f),  and  (6h)  we  have  for  the  small  aperture  angle  of  arrival 
expressions 


(34) 


2  2.  f 

<a  >  =  ir  L 

0 


[1 


2 

+  -jr  sin  Sjpl  k3  «N(>c)  d< 
K  L  ’ 
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(35) 


Bfi(p) 


2tt2L 


['  +-xsi"Er] 

<  L 


J,(kp)  - 

K  ^(k)  dK 


KP 


(36) 


[1 


+  -y-  sin  ypl  Jq(kp)  K3  §^(k)  d<  -  Bg(p). 


Similar  expressions  were  also  obtained  by  Strobehn.[17]  The  corre¬ 
sponding  equations  for  the  spherical  wave  case  can  be  found  by  making 
the  substitution 


(37) 


1  +  -j-  sin 
k  L 


A 


L  r 


Vk) 


k  cos  (ilfcsfc-) 


in  Eqs.  ( 34 ) - ( 36 ) .  These  expressions  are  suitable  for  numeric 
evaluation  since  ^(k)  approaches  zero  rapidly  for  a  finite  value 
of  k  so  that  the  upper  limit  of  the  integrals  can  be  made  finite. 


As  noted  in  Appendix  B,  when  the  outer  scale,  L0  becomes 
infinite,  Eq.  (33)  reduces  to  the  Kolmogorov  spectrum  with  inner 
scale,  and  Eqs.  (32)  and  (34)-(36)  can  be  integrated  ana¬ 

lytically.  Thus,  defining 


D  = 


2  2 

p 

and  g  = 


we  have 


2b 


(38) 


Bg(P)  *  1.028  L  CN2  <J/3 


6  lF1^6,2’"g) 


V. 


(40)  <a2>  *  <e2>  =  L  CN2  Km1/3  (.904  +  1.088  D'1/6  cos 


For  large  p,  the  asymptotic  forms  are 


W)  8B(p)  4  2.43  L  CN2  p’173 

(42)  Bo(0)  4  1.62  L  CN2  o‘1/3. 


The  small  aperture  angle  of  arrival  mean  square  and  correlation 
function  were  calculated  for  many  cases.  These  will  now  be  presented. 
Equations  (38),  (39),  and  (40)  were  evaluated  for  several  typical 
values  of  range  and  inner  scale.  The  results  are  shown  in  Figs.  7-34. 
The  parameters  used  were  chosen  to  agree  with  those  that  might  be 
present  at  the  RADC  Laser  Propagation  Range.  For  example  the  ranges  of 
304.8  m,  914.4  m  and  1524  m  correspond  respectively  to  1000  ft,  3000  ft 
and  5000  ft,  all  ranges  used  at  the  RADC  Range. 
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First  the  mean  square  angle  of  arrival  is  considered.  Figures 
7  and  8  show  it  as  a  function  of  range,  first  for  two  inner  scales  and 
infinite  outer  scale  as  given  in  Eq.  (40),  and  second  for  three  different 
outer  scales  taking  the  inner  scale  as  1  rrm  and  numerically  integrating 
Eqs.  (33)  and  (34).  For  simplicity  the  mean  square  is  normalized  to 
the  refractive  index  structure  parameter,  C^.  it  is  noted  that  the 
mean  square  increases  as  the  outer  scale  increases  and  as  the  inner 
scale  decreases .  In  either  case  the  range  of  turbulent  spatial  fre¬ 
quencies  interacting  with  the  light  beam  has  increased,  thus  increasing 
the  effect.  Of  the  two  scales  the  inner  scale  has  the  larger  effect  on 
the  arrival  angle  flucuations,  because,  while  the  smaller  scale  eddies 
are  less  prevalent,  they  scatter  light  at  the  greatest  angles. 

Typical  angle  of  arrival  correlation  curves  are  shown  in  Figs. 

9,  10,  and  11,  evaluated  from  Eqs.  (34),  (35),  and  (36).  They  show 
the  effect  of  inner  scale  on  the  correlation  functions.  Figures  9 
and  10  show  the  autocorrelation  functions  Ba(p)/<a^>  and  Bg(p)/<cr> 
for  a  range  of  304.8  m,  infinite  outer  scale  and  inner  scales  of  1  cm 
and  1  mm.  Figure  11  shows  the  left  hand  portion  of  Fig.  10  in  detail 
to  display  the  fine  structure  of  the  initial  roll-offs.  The  rapid 
roll-offs  of  the  autocorrelations,  as  the  separation  is  increased 
slightly  from  zero,  is  evident.  However,  as  the  separation  exceeds 
the  inner  scale,  the  curves  level  off.  The  roll-off  and  break  point 
are  strongly  dependent  on  the  inner  scale. 

The  inner  scale  dependence  can  be  explained.  As  noted  earlier 
for  p  <  ju,  the  phase  structure  function  obeys  a  square  power  law  and 
the  correlations  Bg(p)  and  Ba(p)  are  equal  to  the  mean  square  angle  of 
arrival.  From  a  physical  standpoint  we  note  that  the  inner  scale  is 
a  characteristic  dimension  of  the  smallest  blob  in  a  turbulence  cell. 
Thus,  we  expect  the  angle  of  arrival  correlation  for  p  <  l  to  be  nearly 
the  same  as  the  mean  square  value.  It  is  very  likely,  that  the  two 
regions  considered  on  the  wavefront  were  perturbed  by  the  same  blob 
or  a  closely  correlated  blob.  As  the  separation  is  increased  and 
approaches  the  inner  scale  the  number  of  blobs  increases  rapidly  in 
comparison  to  the  increase  in  separation.  However,  as  the  separation 
becomes  sufficiently  large  the  total  number  of  blobs  masks  the  increase 
in  blobs.  Thus,  the  roll-off  of  the  correlation  slows  down  and  levels 
off. 


The  inner  scale  dependence  can  be  derived  from  the  expressions 
used  earlier.  Since  we  are  dealing  with  an  inner  scale  phenomenon 
the  outer  scale  can  be  ignored  and  the  analytic  expressions  are  used. 
Thus,  for  p  <<  iQ  we  have 


(43) 


BB(p) 

<b2> 


?  ? 

.0206  ►  ‘V 
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Fig.  7.  Small  aperture  mean  square  angle  of  arrival 
versus  range  for  infinite  outer  scale  and 
two  inner  scales. 


i.  01  m 


Small  aperture  elevation  angle  of  arrival  autocorrelation  function 
for  two  inner  scales. 


I 

I 
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Fig.  10.  Small  aperture  azimuth  angle  of  arrival  autocorrelation  function 
for  two  inner  scales. 
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Fig.  11.  Comparison  of  exact  and  approximate  equations  for  the  small 
aperture  azimuth  angle  of  arrival  autocorrelation  function. 


within  a  correction  of  D"  cos  tt/12.  For  p  «  sl0  we  see  that  p 
will  be  less  than  one  and  from  Eq.  (43),  Bg(p)  *  <e2>  as  expected. 

A  comparison  of  the  exact  result  and  the  approximate  result  is  shown 
in  Figs.  11  and  12.  Figure  11  shows  both  cases  of  jl0  =  1  mm  and  1  cm 
and  Fig.  12  is  a  magnified  view  of  the  case  of  ju=  1  cm.  The  structure 
function,  Dm(p)/Cm2,  is  also  shown  in  Fig.  12.  We  see  that  the  approxi¬ 
mation  given  in  Eq.  (43)  breaks  down  around  the  point  where  Dty(p) 
changes  from  a  r2  dependence  to  r2/3  dependence  as  expected. 

The  information  in  Figs.  9,  10,  and  11  is  repeated  in  Figs.  13-16, 
along  with  other  curves  to  show  the  effect  of  varying  the  range.  Instead 
of  plotting  the  curves  on  linear  scale  which  has  a  rapid  variation  near 
zero  separation,  the  urves  are  plotted  on  log-log  paper.  Also,  instead 
of  plotting  the  autocorrelation  which  suppresses  the  range  information, 
the  correlation  function  is  plotted,  normalized  for  convenience  to  the 
turbulent  structure  parameter.  The  fact  that  the  curves  are  identical 
except  for  a  vertical  displacement  indicates  that  the  range  information 
is  primarily  in  the  mean  square  for  infinite  outer  scale. 

The  effects  of  the  outer  scale  obtained  from  computer  evaluation 
of  Eqs.  (33)-(36)  are  shown  in  Figs.  17-20.  The  azimuth  and  elevation 
angle  correlation  are  shown  in  Figs.  17  and  18  for  a  range  of  304.8  m, 
inner  scale  of  1  mm  and  three  values  of  the  outer  scale.  The  information 
is  repeated  in  Figs.  19  and  20  as  log-log  plots  of  the  two  arrival  angle 
correlation  functions  normalized  to  the  turbulence  structure  parameter, 


The  effect  of  the  outer  scale  is  shown  quite  dramatically.  The 
elevation  angle  of  arrival,  Ba(p)/<az>  in  Fig.  17  shows  a  slightly  negative 
correlation  for  p  ■  L0.  This  negative  correlation  is  sufficiently  large 
that  it  is  not  believed  to  be  an  artifact  of  the  computation.  Further 
showing  the  effect  of  outer  scale,  the  azimuth  angle  of  arrival  correlation 
Bg(p)/<or>  shown  in  Fig.  18  falls  below  5%  of  the  mean  value  for  p  ■  L0. 

The  effects  of  range  with  a  finite  outer  scale  are  shown  in 
Figs.  21  and  22  where  we  see  the  normalized  correlation  functions 
plotted  for  an  outer  scale  of  1  m,  and  inner  scale  of  1  mm  and  ranges 
of  304.8  m  ft,  914.4  m  ft  and  1524  m  ft.  Figures  23  and  24  contain  the 
same  information  except  with  an  inner  scale  of  1  cm.  It  is  noted  that 
In  these  cases  the  outer  scale  dependence  is  almost  entirely  determined 
by  the  variation  of  the  mean  square  arrival  angle  caused  by  the  outer 
scale.  Thus,  the  correlation  function  is  almost  directly  proportional 
to  range. 

The  proceeding  graphs  have  all  been  for  a  plane  wave.  The 
equivalent  graphs  are  also  presented  for  a  spherical  wave  In  Figs. 

25-34.  Thus  Figs.  25  and  26  give  the  mean  square  arrival  angle  for 
two  outer  scales  for  1  mm  and  1  cm  inner  scales  respectively.  Figures 
27  and  28  give  the  normalized  arrival  angle  correlation  functions  for 
three  ranges  for  an  infinite  outer  scale  and  inner  scale  of  1  mm. 
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Fig.  13.  Small  aperture  elevation  angle  of  arrival  correlation  versus 
separation  for  plane  wave  input,  L  =  »,  t  =  0.01  m. 


5.  Small  aperture  elevation  angle  of  arrival  correlation  versus 
separation  for  plane  wave  input,  L  =  a  =  0.001  m. 


1524  m 


Small  aperture  azimuth  angle  of  arrival  correlation  versus  separation 
and  large  aperture  mean  square  angle  of  arrival  versus  aperture  radius 
for  plane  wave  input,  L  =00,0=  q.001  m. 
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Fig.  17.  Small  aperture  elevation  angle  of  arrival  auto¬ 
correlation  versus  separation  for  plane  wave 
input  and  three  outer  scales. 


Fig.  18.  Small  aperture  azimuth  angle  of  arrival  auto¬ 
correlation  versus  separation  for  plane  wave 
input  and  three  outer  scales. 
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Small  aperture  elevation  angle  of  arrival  correlation  versus  separation 
for  plane  wave  input  and  three  outer  scales. 


Small  aperture  azimuth  angle  of  arrival  correlation  versus  separation 
and  large  aperture  mean  square  angle  of  arrival  versus  aperture  radius 
for  plane  wave  input  and  three  outer  scales. 
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Fig.  21.  Small  aperture  elevation  angle  of  arrival  correlation  versus  separation  for 
plane  wave  input  L  =  1  m,  a  =  .001  m. 
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Fig.  23.  Small  aperture  elevation  angle  of  arrival  correlation  versus  separation  for 
plane  wave  input,  L0  =  1  m,  £n  =  0.01  m. 


Small  aperture  azimuth  angle  of  arrival  correlation  versus  separation  and  large  aperture  mean 
square  angle  of  arrival  versus  aperture  radius  for  plane  wave  input,  L0  =  1  m,  i0  =  0.01  m. 
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Fig.  25.  Small  aperture  mean  square  angle  of  arrival 
for  spherical  wave  input,  finite  outer  scale 
and  two  inner  scales. 
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Fig.  26.  Small  aperture  mean  square  angle  of  arrival 
for  spherical  wave  input,  infinite  outer 
scale,  and  two  inner  scales. 
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Fig.  27.  Small  aperture  elevation  angle  of  arrival  correlation  versus 
separation  for  spherical  wave  input,  L  =  »,  i  =  0.001  m. 
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Fig.  28.  Small  aperture  azimuth  angle  of  arrival  correlation  and  large  aperture  mean  square 
angle  of  arrival  versus  aoerture  radius  for  spherical  wave  input,  L  =  °°*  £0  = 


Figures  29  and  30  give  the  same  data  except  that  the  inner  scale  is 
1  cm.  Figures  31  and  32  give  the  normalized  arrival  angle  correlations 
for  a  1  m  outer  scale  and  1  mm  inner  scale  for  the  three  different 
ranges  and  Figs.  33  and  34  hold  for  1  m  outer  scale,  1  cm  inner  scale 
and  the  three  different  ranges.  These  were  all  computed  the  same  as  the 
plane  wave  data,  except  with  Eq.  (32)  replacing  (31). 

The  accuracy  of  the  computer  evaluations  of  Eqs .  (33)- (36)  was 
checked  by  evaluating  them  for  the  case  where  L0  =  »  and  comparing  with 
the  closed  form  analytic  expressions  in  Eqs.  (38)-(40).  Several  orders 
of  Gaussian  quadrature  integration  were  tried  and  compared  to  the 
closed  form  results  as  shown  in  Fig.  35.  We  see  that  the  96-point 
Gaussian  quadrature  integration  used  in  the  previous  calculations  provides 
good  agreement. 

Extensive  results  are  also  available  for  the  large  aperture  arrival 
angle  variance.  Figures  36  and  37  are  plots  of  mean  square  arrival  angle 
normalized  to  turbulence  structure  parameter,  plotted  as  a  function  of 
aperture  radius.  The  numbers  were  obtained  by  numerically  integrating 
Eq.  (19)  with  the  expression  in  Eq.  (33)  used  for  <?n(k).  In  Fig.  36 
curves  are  shown  for  plane  waves  for  an  inner  scale  of  1  inn,  a  range  of 
304.8  m  ft  and  three  values  of  outer  scale.  Figure  37  shows  the  effect 
of  range  for  an  outer  scale  of  1  m  and  an  inner  scale  of  1  mm.  These 
curves  show  changes  with  outer  scale  and  range  similar  to  the  small 
aperture  arrival  angle  correlation  functions. 

Verifying  the  accuracy  for  these  results  is  much  more  difficult 
and  time  consuming  since  four  numeric  integrations  are  used.  However, 
increasing  the  order  of  integrations  used  yielded  results  which  differed 
in  the  third  and  higher  order  digits.  Thus,  these  results  are  believed 
to  be  accurate  to  3  digits.  Higher  accuracy  can  be  obtained  by  in¬ 
creasing  the  Gaussian  integration  orders. 

Further  data  is  also  available  on  large  aperture  mean  square 
arrival  angle,  because,  to  a  very  good  degree  of  approximation,  the 
large  aperture  mean  square  arrival  angle  taken  as  a  function  of  aperture 
radius,  is  equal  to  the  small  aperture  azimuth  arrival  angle  correlation 
function,  Bg(p)  taken  as  a  function  of  point  separation,  p.  This 
synonymity  is  discussed  in  detail  in  the  next  section.  However  the  results 
are  used  here  to  indicate  the  further  data  on  large  aperture  arrival 
angle  mean  square  available.  Thus  the  even  numbered  figures  between 
Fig.  14  and  Fig.  24  can  all  be  retitled  to  give  large  aperture  mean 
square  angle  of  arrival  versus  radius  for  the  particular  case  for 
plane  waves,  and  the  even  numbered  figures  between  Fig.  27  and  Fig. 

34  can  all  be  retitled  to  give  the  mean  square  arrival  angle  versus 
radius  for  the  particular  cases  for  spherical  waves.  Thus  data  is 
available  for  a  large  number  of  cases  for  the  large  aperture  situation 
also.  The  discussion  of  the  curves  parallels  that  for  the  small  aperture 
case . 
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Fig.  29.  Small  aperture  elevation  angle  of  arrival  correlation  versus  separation  for 
spherical  wave  input,  L  =  ®,  z  =  0.01  m. 
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Small  aperture  elevation  angle  of  arrival  correlation  versus  separation 
for  spherical  wave  input,  L  =  1  m,  fc  =  0.001  m. 
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Fig.  33.  Small  aperture  elevation  angle  of  arrival  correlation  versus  separation 
for  spherical  wave  input,  L  =  1  m,  l  =  0.01  m. 
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of  arrival  correlation  functions. 


This  concludes  the  presentation  of  data.  Large  and  small  aperture 
angle  of  arrival  data  predicting  experimental  results  have  been 
calculated  for  a  range  of  typical  cases.  The  effects  of  range,  inner 
and  outer  scales,  separation  and  aperture  size  have  all  been  demonstrated. 

Since  all  of  the  equations  derived  in  this  report  are  functions 
of  Dc(p),  the  validity  of  these  results  depends  on  the  validity  of 
the  Rytov  approximation.  This  topic  will  now  be  briefly  considered. 

The  question  has  been  the  subject  of  many  discussions[18-26]  and  the 
criterion  which  seems  appropriate  at  this  time  to  the  calculations  in 
this  report  is  Ds(p)  «  ir.  This  result  was  reached  by  Tatarski  from 
calculations  of  the  second  order  term  of  Ds(p).[18]  A  similar  result 
was  also  reached  by  DeWo'if.[19]  To  quantify  this  result  we  consider 
p  much  greater  than  i0  and  set  [5,(47.37)] 


(44)  Ds(p)  =  2.92  CN2  k2  p5/3  L  «  ir. 


For  convenience,  replace  the  "much  less  than  tt"  with  ir/10.  Thus,  we 
can  write  Eq.  (44)  as 


(45) 


p 


This  is  shown  plotted  in  Fig.  38. 

In  Fig.  38  several  range,  Cfj2L,  combinations  are  shown.  We  see 
that  for  the  longest  range  used  earlier,  1524  m  and  for  strong 
turbulence  C^2  =  1 0“ 1 ^  m-2'3,  the  maximum  permissible  separation,  p, 
is  approximately  9.2  cm.  This  means  that  the  angle  of  arrival  two-point 
correlation  is  valid  for  p  <  9.2  cm  and  the  large  aperture  results  are 
valid  for  R  <  9.2  cm.  As  the  range  and  turbulence  levels  decrease  we 
note  that  the  restrictions  are  relaxed. 

It  should  also  be  noted  that  the  form  used  for  Ds(p)  does  not 
include  saturation  effects  due  to  the  outer  scale.  Since  the  validity 
criterion  was  derived  independently  of  the  exact  form  of  the  refractive 
index  fluctuations  saturation  effects  may  extend  the  region  of  validity. 
Furthermore,  it  has  been  suggested  the  region  of  validity  may  be  larger 
than  calculated. [5], [27]  Measurements  by  Bouricius  and  Clifford  at 
.632y[27]  suggest  this  result  and  indicate  that  more  study  is  required. 
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Fig.  38.  Areas  of  validity  for  Rytov  approximation. 


To  summarize  the  work  presented  in  this  section,  data  were 
presented  predicting  results  of  various  angle  of  arrival  experiments 
for  various  conditions.  First  the  formulae  to  be  used  for  the  phase 
structure  function  and  refractive  index  turbulence  spectra  were  pre¬ 
sented.  Then  expressions  given  in  previous  sections  tor  small  aperture 
arrival  correlation  functions  were  evaluated,  and  plotted.  For 
infinite  outer  scale  the  expressions  could  be  written  as  closed  form 
expressions,  for  finite  outer  scale  they  had  to  be  evaluated  numerically. 
Both  elevation  and  azimuth  angle  correlation  functions,  respectively 
B  (p)  and  Bg(p),  were  calculated  and  graphed  for  plane  waves  for  ranges 
304.8  m,  914.4  m  and  1524  m,  outer  scales  of  infinity,  1  m  and  0.159  m, 
and  inner  scales  of  0.001  m  and  0.01  m.  This  data  was  then  repeated 
for  spherical  waves.  The  large  aperture  mean  square  arrival  angle  was 
also  evaluated  from  a  previously  presented  formula.  Only  a  few  cases 
were  evaluated  because  of  a  synonymity  which  allowed  the  small  aperture 
azimuth  arrival  angle  correlation  curves  as  a  function  of  separation,  p, 
to  be  relabled  with  large  aperture  mean  square  arrival  angle  as  a 
function  of  aperture  radius,  R,  thus  allowing  many  of  the  results  to 
serve  double  duty.  The  section  was  concluded  with  a  discussion  of  the 
Rytov  approximation  on  which  the  phase  structure  function  was  based. 


VI.  ANGLE  OF  ARRIVAL  SYNONYMITY 

This  section  deals  explicitly  with  a  relationship  between  the 
angle  of  arrival  variance  for  large  aperture  receiver  as  a  function 
of  aperture  radius,  R,  and  the  small  aperture  azimuth  correlation 
function,  Bg(p)  taken  as  a  function  of  p.  The  particular  relationship 
is  that  they  are  synonymous  over  a  sizable  range  of  values  of  the  argument. 
These  two  situations  are  shown  in  Fig.  39  to  facilitate  the  comparison 
between  the  two.  Figure  39a  shows  a  standard  large  aperture  arrival 
angle  measurement  situation.  An  atmospherically  degraded  beam  from 
a  point  source  is  focussed  to  a  degraded  spot.  The  information  about 
angle  of  arrival  is  obtained  from  the  transverse  motion  of  the  spot, 
the  arrival  angle  being  proportional  to  the  spot  displacement. 

Figure  39b  shows  a  possible  situation  for  the  measurement  of  small 
aperture  arrival  angle  correlation  function.  The  angles  of  arrival 
at  two  pinholes  are  determined  independently  by  monitoring  the  two 
pinhole  diffraction  patterns,  the  motion  of  the  center  of  each  pattern 
being  proportional  to  the  individual  arrival  angle.  In  the  case  at 
hand,  only  the  x  component  of  motion  is  of  interest.  The  pinholes  are 
sufficiently  small  that  the  amplitude  is  constant  over  the  hole.  Then 
the  quantity  of  interest  is  the  correlation  of  arrival  angles  between 
these  two  pinholes,  the  correlation  being  a  function  of  p,  the  separation 
of  the  pinholes.  In  the  case  of  interest  the  pinholes  are  separated  in 
the  direction  perpendicular  to  the  direction  in  which  the  spot  motion 
is  being  monitored,  (the  x  direction  in  Fig.  39). 
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Fig.  39.  Comparison  of  the  definitions  of  the  large  aperture 
mean  square  angle  of  arrival  and  the  small  aperture 
angle  of  arrival  correlation. 
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Comparing  Eqs.  (46)  and  (47)  we  see  that  both  expressions  can  be 
written  in  the  form 


(48)  4it^L  /  T(<)  F ( tc )  4>^(k)  die 


where  T(«)  is  the  filter  function  peculiar  to  the  particular  operation, 
F(k)  identifies  the  type  of  wave,  and  ^(k)  is  the  refractive  index 
turbulence  spectrum.  Denoting  the  particular  small  and  large  aperture 
filter  functions  by  SAUp)  and  LA(kR)  respectively,  we  have 


(49)  SA(<p)  =  J-|  ( <p  )/<p 

R  R  ir  _ 

(50)  LA(kR)  =  -  ( 1  6/ttR®<^)  |  r^  dr  |  p^  dp  j  cos  y  [l-J0(tcir^+p^-2rpCOS  y)]dy 

ooo 


Equations  (49)  and  (50)  describe  the  manner  in  which  the  receiving 
aperture  interprets  the  angle  of  arrival.  Graphs  of  these  functions 
are  shown  in  Fig.  41.  For  convenience  in  comparison  we  shall  take 
them  both  to  be  functions  of  the  same  variable  to  be  called  R.  We 
see  in  this  figure  that  the  two  filter  functions  are  synonymous  over 
a  wide  range  of  <R. 

The  synonymity  is  limited  with  respect  to  the  range  of  the  variable 
R  over  which  it  applies.  This  occurs  because  of  the  relative  size 
of  the  filter  function  and  the  other  functions,  F(k),$n(k)  k3  in  the 
frequency  integral.  If  the  other  functions  are  of  significant  size 
only  when  the  filter  functions  are  almost  identical,  i.e.,  in  the 
region,  0  <kR<4,  then  the  identity  will  hold.  If  the  other  functions 
have  large  contributions  for  <R>4  then  the  synonymity  will  break  down. 

To  more  quantitatively  define  the  regions  of  applicability  for  the 
angle  of  arrival  synonymity  we  define  e(R),  the  error  between  the  two 
functions  relative  to  the  small  aperture  correlation  function, 


,  ,  MR)  -  <ax S 

(51)  e(R)  =  BQ(R) 
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arrival  filter  functions. 


The  numerator  of  Eq.  (51)  is  given  bv  an  integral,  I d( kR)  of  The 
difference  between  Eqs.  (46)  and  (47).  Using  the  difference  function 
a(<R)  given  in  Eq.  (52)  and  shown  in  Fig.  42, 


(52)  a(kR)  =  SA(<R)  -  LA(kR) 


the  numerator  integral  is 


(53) 


Iq(kR)  =  4tt^L 


00 

a(kR)  Fs(k)  <I>hj(k) 

0 


3 

K 


die. 


For  a  plane  wave,  Fs(<)  varies  between  1  and  2  and  has  little  effect, 
thus,  making  a(kR)  and  «3  4>N(<)  the  two  functions  of  most  interest. 

We  now  proceed  to  examine  e(R)  qualitatively,  but  in  some  detail 
to  determine  the  approximate  functional  dependence.  This  is  to  be 
done  by  examining  the  curves  for  SA(kR),  LA(kR),  a(kR),  and  <3  $»(K) 
and  the  various  ways  they  can  combine  for  various  values  of  R.  These 
curves  are  plotted  quite  precisely  in  Figs.  41,  42,  and  43,  using  for 
$n(k)  the  Von  Karman  spectrum  given  in  Eq.  (33).  The  various  functions 
are  compared  in  Fig.  44  for  several  values  of  R.  The  effect  of 
choosing  various  values  of  R  is  to  select  particular  horizontal  scales 
for  the  filter  functions  SA(kR),  LA(kR)  and  a(kR).  Thus,  as  seen 
in  Fig.  41,  LA(kR)  goes  essentially  to  zero  at  <  =  4.4/R  and  SA(kR) 
has  its  first  zero  at  k  =  3.9/R.  Fig.  42  shows  that  the  difference 
function,  a(kR)  has  its  first  zero  at  <  =  3.4/R,  pretty  much  in  the 
same  region.  Examining  curves  in  Fig.  44a  and  we  see  that  for 
R  «  £0  The  small  and  large  aperture  filter  functions  are  nearly 
constant  over  the  complete  range  of  k3  4^(k)  so  that  Bg[R)  and  <avZ(R)> 
are  identical,  each  being  given  by  that  integral  over  <3  $«(<).  For 
R  =  4.4/s,0,  comparing  Figs.  44a  and  44c,  we  see  that  both  LA(<R)  and 
SA(kR)  go  to  zero  at  about  the  same  place  that  k3  4>m(k)  becomes  very 
small.  For  that  case  both  curves  are  finite  over  the  range  0-l/Lo, 
the  product  at  any  point  being  less  than  the  value  of  either  curve. 

The  integral  over  the  product  is  then  starting  to  decrease.  This  is 
the  region  of  R  where  the  curves  of  Bg(p)  and  <ax2(R)>  start  to 
decrease,  as  shown  on  graphs  of  these  functions.  Further  the  value 
of  the  error  f  a(kR)  k3  «>n(k)  d<  has  been  increasing  with  R 
and  has  possibly  a  local  maximum  because  a(kR)  has  a  negative  portion 
where  k3  ^(k)  is  still  finite. 
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Fig.  42.  Large  and  small  aperture  angle  of  arrival  difference 
function  versus  spatial  frequency. 


10" 


fig.  44.  Comparison  of  large  and  small  aperture  filter 
functions  with  the  modified  refractive  index 
fluctuation  spectrum. 
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For  R  in  the  region  of  L0,  comparing  Figs.  44a  and  44d  shows 
that  a(kR)  and  k3  $n(k)  both  have  large  positive  portions  in  the 
region  of  k  =  1/L0.  One  would  expect  a  local  maximum  in  the  error 
around  this  region.  As  R  increases  slightly  past  L0,  the  first  zero 
of  the  a(kR)  curve  approximately  coincides  with  the  peak  of  the 
k3  <&n(k)  curve  so  that  the  product  is  positive  for  k<1/L0  and  negative 
for  <  >  l/L0.  One  would  expect  the  error  to  decrease  for  this  region  of 
R.  Finally  as  R  becomes  still  larger,  the  error  curve  will  oscillate. 
This  expected  qualitative  behavior  is  sketched  in  Fig.  45. 


Fig.  45.  Integrated  large  and  small  aperture  angle 
of  arrival  difference  function. 


One  can  further  examine  the  relative  error,  e,  defined  in  Eq.  (51). 
In  particular  for  the  case  where  R  <  L0  the  product  of  LA(kR)  and 
k3  4>n(k)  becomes  quite  small,  being  limited  for  small  k  by  k3  *m(k)  and  for 
large  k  by  LA(kR).  Then  <ax2(R)>  in  Eq.  (51)  is  relatively  small  and 
e(R)  approaches  unity.  Thus  the  assymptotic  value  of |e ( R) J i s  1.00, 
i.e.,  100%  error.  The  expected  qualitative  behavior  for  e(R)  is  then 
plotted  using  the  data  on  a(<R),  SA(kR)  and  the  assymptotic  value.  The 
plot  is  given  in  Fig.  46. 
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Fig.  46.  Angle  of  arrival  synonymity  error  function 
showing  asymptotic  behavior. 


These  qualitative  results  are  confirmed  by  the  quantitative 
numeric  evaluation  in  Eq.  (53),  (including  Fs(<))  shown  in  Fig.  47. 
These  numeric  calculations  required  considerable  computation  time  so 
only  enough  points  were  calculated  to  show  the  nature  of  the  curves. 
Thus,  some  oscillations  have  more  than  likely  been  suppressed  and  the 
final  value  of  e(R)  was  not  reached.  Accuracy  in  these  calculations 
was  difficult  to  obtain,  however,  the  results  ar''  certainly  correct 
to  at  least  order  of  magnitude  and  demonstrate  tne  behavior  predicted 
earlier. 

Three  cases  are  shown  in  Fig.  47  corresponding  to  earlier 
calculations.  The  first  characteristic  we  note  is  that  e(R)  is 
fairly  small  for  R  <  L0.  This  implies  that  for  aperture  diameters 
less  than  twice  the  outer  scale,  the  large  aperture  mean  square 
angle  of  arrival  can  be  approximated  by  the  small  aperture  azimuth 
correlation  function.  From  Fig.  47  we  see  that  this  error  will  be 
on  the  order  of  10%.  Second,  we  note  the  range  has  very  little 
effect  as  is  to  be  expected.  What  effect  the  range,  L,  does  have  is 
incorporated  in  Fs(<). 

To  summarize,  in  this  section,  we  have  stated  that  the  large 
aperture  angle  of  arrival  mean  rquare  value  and  the  small  aperture 
azimuthal  angle  of  arrival  correlation  function  are  synonymous.  The 
functional  dependence  was  demonstrated  to  be  very  nearly  equal  for 
typical  cases.  The  relationship  was  examined  in  the  spectral  domain 
both  qualitatively  and  quantitatively  showing  that  the  results  should 
be  expected  to  be  well  within  ten  percent  for  values  of  the  argument 
up  to  the  outer  scale. 
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Fig.  47.  Angle  of  arrival  synonymity  error  function  for  two 
outer  scales  and  two  ranges. 
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VII.  SUMMARY 


In  this  report  the  subject  of  angle  of  arrival  of  light  beams 
degraded  by  propagation  through  atmospheric  turbulence  has  been 
considered.  Data  on  arrival  angle  experiments  was  calculated  for  two 
types  of  experiments  for  a  wide  variety  of  cases  and  an  equivalence 
between  two  types  of  angle  of  arrival  experiments  was  demonstrated.  In 
Chapter  II  the  subject  was  Introduced  and  reviewed,  and  the  literature 
enumerated.  It  was  pointed  out  that  there  were  two  approaches  to 
angle  of  arrival  considerations,  one  using  a  pair  of  pinhole  apertures 
and  the  other  using  a  single  large  aperture.  In  Chapter  III  the 
area  of  small  aperture  angle  of  arrival  calculations  was  reviewed  In 
detail,  the  basic  theoretical  expressions  being  derived.  The  calculation 
of  large  aperture  angle  of  arrival  was  discussed  In  Chapter  IV.  In 
Chapter  V  the  analytic  expressions  derived  in  Chapters  III  and  IV  were 
evaluated  using  a  particular  expression  for  the  refractive  Index  spatial 
spectrum,  namely  the  Von  Karman  spectrum.  The  results  were  presented 
graphically  for  the  two  angle  of  arrival  situations  for  variety  of  cases 
Including  various  Inner  and  outer  scales.  The  Itytov  approximation  used 
In  the  derivation  of  the  phase  structure  function  used  was  considered 
and  a  graph  given  showing  the  expected  range  of  validity  the  calculations. 
In  Chapter  VI  a  synonymity  was  presented  between  the  large  aperture 
angle  of  arrival  mean  square  value  expressed  as  a  function  of  aperture 
radius  and  the  small  aperture  azimuth  angle  of  arrival  correlation 
function  taken  as  a  function  of  aperture  separation.  This  synonymity 
v;es  used  In  Chapter  V  to  simplify  calculations  and  also  offers  the 
possibility  of  simplification  of  experimental  measurements.  The 
limitations  on  the  synonymity  were  demonstrated. 
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APPENDIX  A 

DISCUSSION  OF  REFRACTIVE  INDEX  FLUCTUATION  SPECTRA 


Various  forms  of  the  refractive  index  fluctuation  spectra 
have  been  used  to  find  solutions  of  the  wave  equation  for  an  inhomo¬ 
geneous  random  media.  Two  forms  which  have  been  used  recently  are 
compared  in  this  appendix.  The  choice  of  one  of  the  constants  is 
explained  and  several  interesting  characteristics  noted.  The 
derivations  of  these  spectra  are  not  presented  since  they  are  outside 
the  scope  of  this  report. 

Kolmogorov  theory  predicts  that  the  refractive  index  structure 
functions  should  have  the  following  asymptotic  behavior 


0  <  r  «  jl 

o 

*0  <,:  r  <<r  Lo 


where 


(Al) 


Wr> 


2  .2/3  ,r  x2 


<£  l, 
N  o 


r—y 

V 


C2  r2/3 

ln  r 
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CN  =  structure  constant  of  the  turbulence. 


It  can  be  shown  that  the  Kolmogorov  spectrum  as  used  by  Tatarski  in  [5] 
and  given  in  Eq.  (A2) 
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has  a  structure  function  of  the  form 
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where  -| F-j  signifies  a  confluent  hypergeometric  function.  Further, 

It  can  be  shown  the  confluent  hypergeometric  function  has  the  proper 
asymptotic  behavior  such  that  Eq.  (Al)  is  satisfied. 

Another  form  of  $«(«)  is  the  Von  Karmann  spectrum  which  includes 
outer  and  Inner  scale  effects.  Several  algebraic  forms  have  been  used 
but  they  can  be  generalized  to 


(A4) 

wi  th 


$n(k)  = 


2.  2 

o  -  K  /  K 

.033  Cjj  e  m 

— rTT^TIW 

(k  +  <0) 


where  A  is  a  normalization  constant.  The  value  of  A  has  been  chosen 

as  1  in  [17]  and  2ir  In  [28].  For  L  =  «,  Eq.  (A4)  reduces  to  Eq.  (A3)  as 

expected.  0 

The  value  of  A  used  in  this  report  is  1.  The  reason  for  this 
choice  will  now  be  explained.  We  note  that  when  A  *  2ir,  the  effect  Is 

to  normalize  L  to  a  new  value  L'  such  that 


(A5) 


The  physical  interpretation  of  Lq  is  that  It  is  a  characteristic 
dimension  of  saturation  effects.  To  examine  this  situation  the  refractive 
Index  spectrum  and  structure  function  calculated  from  Eq.  (A3)  and  Eq.  (A4) 
and  the  resulting  phase  structure  function  are  shown  in  Figs.  A1-A3  re¬ 
spectively.  Dashed  lines  have  been  added  at  p  =  l/«  =  L  for  =  1 
and  =  1/2tt.  0  0 

Two  Interesting  characteristics  are  noted.  First,  all  of 
the  curves  tend  to  saturate  at  a  value  predicted  approximately  by 
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Fig.  A1 .  Refrartlve  Index  fluctuation  spectrum  versus 
spatial  frequency. 
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the  Kolmogorov  spectrum  at  p  =  l/«  =  L0.  Secondly,  the  break  point  of 
the  saturation  occurs  at  the  same  point.  These  two  characteristics 
are  consistent  with  the  original  definition  of  the  outer  scale,  L0, 
and  the  use  of  A  =  2ir  merely  clouds  the  issue.  Thus,  for  this  report 
A  is  chosen  to  be  one. 

Another  interesting  characteristic  is  seen  in  Fig.  A3.  The 
slope  of  the  phase  structure  function  in  the  inertial  subrange  is  not 
unique.  Since  the  effect  of  the  outer  scale  is  to  cause  saturation 
the  slope  must  change  continuously  from  that  of  square  law  to  zero. 
Only  for  the  special  case  of  L0  =  «  is  the  slope  a  constant.  There¬ 
fore,  the  outer  scale  makes  it  difficult  to  predict  a  unique  slope  of 
the  phase  structure  function  in  the  inertial  subrange. 

Concluding,  two  refractive  index  fluctuation  spectra  have  been 
described  and  the  criteria  used  to  determine  one  of  the  constants 
explained.  Several  interesting  characteristics  were  demonstrated. 


APPENDIX  B 

COMPUTER  PROGRAM  USED  FOR  NUMERIC  CALCULATIONS 


The  computer  programs  used  for  numeric  calculations  in  this 
report  are  presented  in  this  appendix.  Each  program  will  be  described 
briefly  and  data  input  explained.  A  listing  of  the  source  deck  is 
also  presented.  All  programs  were  written  for  use  on  an  IBM  360/75 
computer  in  FORTRAN  IV(G).  Some  subroutines  are  contained  in  the  IBM 
Scientific  Subroutine  Package  (360-CM-03X)  Version  III. 


1 .  Small  Aperture  Angle  of  Arrival  Calculations 

Three  programs  were  written  to  make  the  small  aperture  calcu¬ 
lations.  The  first  program  evaluates  the  analytic  functions  for  the 
angle  of  arrival  mean  square  and  2  pt.  correlation  function,  using  the 
Kolmogorov  spectrum  and  plane  wave  input.  The  second  program 
evaluates  the  integral  expressions  using  the  Von  Karmann  spectrum  and 
plane  wave  input.  The  third  program  evaluates  the  integral  expressions 
for  spherical  wave  input.  These  programs  also  calculate  the  refractive 
index  fluctuation  spectrum  and  structure  function,  the  phase  structure 
function  and  the  wave  structure  function.  Each  program  has  numerous 
comment  cards  to  describe  the  program  flow  and  the  variable  names 
correspond  closely  to  variables  used  in  the  report. 

(a).  The  computer  program  used  to  calculate  the  angle  of 
arrival  correlation  functions  for  the  Kolmogorov  spectrum  is  shown  in 
listing  (Bl).  It  evaluates  the  confluent  hypergeometric  functions  for 
pure  real  or  pure  imaginary  arguments  using  formulae  taken  from  Tables 
of  Functions,  Jahnke  and  Emde. 

The  first  card  in  the  input  data  specifies  the  parameters,  t0, 

L,  and  x.  The  variables  are  punched  into  ten  column  fields  in 
decimal  form.  They  are: 


Col umns 

1-10  XLO  =  i  ,  the  inner  scale  of  the  turbulence  in  meters 

11-20  XL  =  L,  the  range  in  meters 

21-30  XLAMB  =  x,  the  wavelength  in  meters. 


If  XL  =  0,  execution  is  terminated.  The  following  card  specifies 
the  starting  point,  increment,  and  number  of  increments  for  the  separ¬ 
ation  distance.  These  variables  are: 


84 


1  LISTING  B-l 

I 

C  10.6  PROPAGATION  PHENOMENON  PROGRAM  **** 

1C  *****  SMALL  APERTURE  ANGLE  OF  ARRIVAL  CALCULATION  FOR 
C  THE  KOLMOGOROV  SPECTRUM  ***** 

C 

I  COMPLEX  G»GO 

REAL*8  R0YP2*XDYP2*RDYP3*XDYP3*  81 »B2 
C  INITIALIZE  NORMALIZATION  CONSTANTS 
BAN0KM*1. 

IBBNORM-1. 

BAG0RM=1. 

BBGORM-1. 

fC  READ  DATA 

25  READ! 5 »9C0)XL0»XL»XLAMB 
IF(XLO.EQ.O. ICALL  EXIT 

150  REAL) (5  *900 )  RHOS»  DELRH* XNR 

IFUNR.EQ.O. )G0  TO  25 
900  FORMAT (3EI0. 3) 

WRI TE! 6«1001 ) 

(  1001  FORMAT  I *1*1 

(  C  CALCULATE  CONSTANTS 

XM6.2821/XLAMB 
>  XKM=5.  94/XLO 

D«XKM**2*XL/XK 

NR*XNR+1 

A1«1.088*XL*X KM**. 333333 
A2«D**!-. 1666667) 

I  A3«D**i-l. 166667) 

BI*DCOS ( 3* 1415927/12 • I 
B2*DSIN(3. 1415927/12.) 

WRI TEI 6» 1002 ) XKvXKMt U 

1002  FORMAT!/*  XK*  *« lPE10.3t5Xt *XKM>  • , E10.3.5X, *D»  *.E10.3) 
C  BEGIN  ITERATION  ON  KHO 
DO  100  1*1 tNR 
RHO*RHOS+DELRH*FLUAT ( I-l) 

G«XKM**2*RH0**2/4. 

XGD* ( RH0**2*XK ) / ( 4 .*XL ) 

GD* (0. *  1. )*XGD 

C  BRANCH  TO  ASYMPTOTIC  CALCULATION 
IF (XGD.GT. 100. )G0  TO  75 
C 

C  SMALL  RHO  FORMULATION 
C 

C  CALCULATE  CONFLUENT  HYPERGEOMETRIC  FUNCTIONS 
CALL  HYPERG ( . 1666*2. »-G#RHYPl »XHYP 1 ) 

CALL  DHYPER(-.8333.2..GDtRDYP2fXDYP2) 

C  PRINT  VALUES  OF  CONFLUENT  HYPERGEOMETRIC 
WRI TE i 6* 1003IRHYP1 »  XHYPl* RDYP2? XDYP2 
1003  FORMAT !4!E10.3«5X) ) 
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C  CALCULATE  BETA  CORRELATION 

BBETA=Al*(.834*RHYPl*(A2*(  6l*ROYP2+  B2*XOYP2m 
C  CALCULATE  CONFLUENT  HYPERGEOMETRIC  FUNCTIONS 
CALL  HYPERG (1.1666*3. *-G*RHYP2»XHYP2 ) 

CALL  OHYPER I • 16667 » 3 • , GD,RDYP3, XDYP3 ) 

C  PRINT  VALUES  OF  CONFLUENT  HYPERGEOMETRIC 
WRITE  1 6» 1003 )RHYP2fXHYP2tROYP3fXOYP3 
C  CALCULATE  ALPHA  CORRELATION 

BALPHA=BBETA*. 2260*RHO**2*XL*XKM**2. 333* I • 166667*RHYP2+A3*C -B1*XD 
1P3+B2*RDYP3)) 

GO  TO  60 
75  CONTINUE 

ASYMPTOTIC  FORMULATION  FOR  LARGE  RHO 

CALL  HYPERG I. 1666,2. ,-G.RHYPi ,XHYPl ) 

BBETA  = l . 81*XL*XKM**. 3333*RHYPl 
CALL  HYPERGf l • 1666,3 • ,-G,RHYP2, XHYP2 ) 

BALPHA»BBETA  -l.8l*XL*XKM**2. 33333*RH0**2/24.*RHYP2 
60  CONTINUE 

C  NORMALIZE  TO  MEAN  SQUARE  ANGLE  OF  ARRIVAL 
BBN*BBET  A/BBNORM 
BAN*BALPHA/BANORM 
C  PRIN!  RESULTS 

WRITE  1 6, 1000) RHO, BBETA, BALPHA 
WRI TE (6,1004) 

1004  FORMAT (T65, 'CORRELATION  NORMALIZED  TO  MEAN  SQUARE  VALUE  M 
WRITE (6, 1000) RHO, BBN, BAN 
I F ( RHO.NE.O. ) GO  TO  100 
BBNORM»BBETA 
BANORM«B ALPHA 
BBGORMaBBETAG 
BAGORM*BALPHG 
100  CONTINUE 

1000  FORMAT ( T65 , *RHO= ' , 1PE10. 3, 5X, *BBETA*», ElO. 3, 5X,*BALPHA*  *,E10.3) 
GO  TO  50 
END 
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SUBROUTINE  OHYPERt  ALPHA»GAMMA*CZ*RHYP  *IHYP ) 

THIS  SUBROUTINE  CALCULATES  THE  CONFLUENT  HYPERGEQMETRIC  FUNCTION 
ARGUMENTS  ARE  ALPHA, GAMMA,  AND  Cl.  RESULT  IS  (RHYP» IHYP )  (COMPLEX) 

Cl  IS  ASSUMED  TO  BE  PURE  IMAGINARY 

IMPLICIT  COMPLEX (C) . REAL*8  ( A-B,D- I #0-2) 

REAL**  ALPHA, GAMMA, ARGZ*ARZG*  A* B* BMA* S IGN*Z 
C0MPLEX*16  CDZ*CRSUM»CSSUM»CRTERM»CSTERM»CF*DCONJG 
INTEGER  I 

ARG1CZZ )*ATAN2 (AIMAG(CZZ)yREAL(CZZ) ) 

1000  FORMAT  ( *0N0  CONVERGENCE* /•  CHYP**  ,  2UPE10. 3,2X )  ,  5X,  'CHYPO*  •  .2  ( E10. 
13  *2X) ) 

1001  FORMAT  (  *  I*  *  *  I5»5X*  *  ALPHA*  *  *F  10.5, 5X, 'GAMMA*  ' ,  F10. 5  ,5X  *  *CZ*  *«2 
H1PE10.3,2X)*3X,  ’RHYP*  ' *E10.3,5X, • IHYP*  **E10.3) 

ARGZ*ARG(CZ ) 

SMALL  ARGUMENT  EXPANSION 

SET  UPPER  LIMIT  ON  THE  NUMBER  OF  TERMS 
LMT-200 
100  1*2 

CALCULATE  REAL  PART 

RHYP*1. 

IHYP*0. 

Z*AIMAG(CZ ) 

C  TEST  FOR  ZERO 

I  FI  ABS ( Z ) .LT . 1. E-50 ) RETURN 
C  TEST  FOR  LARGE  ARGUMENT 

IF (ABS ( AIMAGtCZ ) ) .GT. 3. )G0  TO  300 

RHYPO*ALPHA/GAMMA*Z*( ALPHA*1 • )/ (GAMMA+1. )*Z/2.*(-l.) 

SRHYP*RHYPO 
120  1*1+2 

C  CHECK  NUMBER  OF  TERMS 

IF ( I.GT.LMT )G0  TO  140 
K*I-1 

RHYP*-RHYPO* ( ALPHA+K- 1 . ) / ( GAMMA +K- 1 • ) *Z/FLOAT (I - 1 ) * ( AL  PHA+K ) / ( GAMM 
1A+K)*Z/FL0AT(I) 

SRHYP*SRHYP+RHYP 
C  TEST  FOR  CONVERGENCE 

IF( DABS ( RHYP-RHYPO) •LT«1*E~05)GO  TO  130 

RHYPO*RHYP 

GO  TO  120 

C  WRITE  ERROR  MESSAGE 
140  NRI TE ( 6* 1000) RHYP.RHYPO *CHYPO 

WRITEt  6*1001 ) I » ALPHA »GAMMA*CZ» RHYP, IHYP 
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130  IHYPU*0 


CALCULATE  IMAGINARY  PART 
RHYP=SRHYP+i. 

1=1  f 

I HYPO=ALPHA/GAMMA*Z 
S IHYP= IHYPO 
220  1=1+2 

CHECK  NUMBER  OF  TERMS 

I F 4 I.GT.LHT1GO  TO  240  > 

IHYP=-IHYPQ*! ALPHA+K-1* 1/tGAMMA+K-l. 1*Z/FL0ATI I-L  )*! ALPHA+K ) / 1 GAMK 
IA+K 1*2/ FLOAT ( 1 1 
SIHYP=1HYP+SIHYP 
TEST  FOR  CONVERGENCE 

IF!UABS!IHYP-IHYP0).LT.I.E-051G0  TO  230 
I  HYPO* I HYP 
GO  TO  220 

WRITE  ERROR  MESSAGE 
240  WRITEtb.iOOO) IHYP, IHYPO, CHYPO 

WRITE (6* 1001 1 1 » ALPHA* GAMMA » CZ »RHYP» I HYP 
230  I HYP=SIHYP 

RETURN  { 

500  CONTINUE 

ASYMPTOTIC  EXPANSION  FOR  LARGE  ARGUMENT 

COZ-CZ 

CRSUM=I. 

CSSUM-1. 

CRTERM=I. 

CSTERM=1. 

C  SET  UPPER  LIMIT  ON  THE  NUMBER  OF  TERMS 
LMT =5 

IFIABSIZJ.GT.  20.)LMT=1 
DO  600  1=1* LMT 
K*  I -L 

CRTERM= ( ALPHA+K 1  ♦ !  AL  PHA-GAMMA+l .+K ) / ( FLOA  T ( I ) *  C-COZ ) ) *CRTERM 

r  p  CIJMsP P  Q  TF QM 

CSTERM* (GAMMA-ALPHA+K ) *  1 i. -ALPHA+K )/ ! FLOAT! I ) *COZ )*CSTERM  2 

CSSUM=CSSUM+CSTERM  [ 

600  CONTINUE 

SIGN»+l.  i 

IF!AIMAG!CZ).LE.0.)SIGN*-1. 

CALL  GAMMBIALPHA* GAMMA* A»B«BMA ) 

ARGZ=3. 141592/2. 

CF»CDABS ICOZ) ♦*( -ALPHA ) *CEXP ( ( 0. *  1 . )* ALPHA* S1GN*ARGZ ) /BMA+CRSUM 
IF! ABS ( Z ) *GT. 20. )G0  TO  650 


CF>CF*  CABS (CZ I I ALPHA-GAMMA )*CEXP( { 0* , I . ) ♦ I  ALPHA-GAMMA ) *ARGZ+C2 ) 
I/A*CSSUM 
650  CONTINUE 
CF»CF*B 

RHYP»« 5* ( OCONJGICF J+CF ) 

IHYP*(0«».5)*(DCONJG(CF )-CF ) 

RETURN 

END 
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SUBROUTINE  HYPERG( ALPHA, GAMMA, CZ,RHYP, IHYP )  i 

THIS  SUBROUTINE  CALCULATES  THE  CONFLUENT  HYPERGEOMETRIC  FUNCTION 
ARGUMENTS  ARE  ALPHA, GAMMA,  AND  CZ.  RESULT  IS  !RHYP»IHYP)  (COMPLEX) 

CZ  IS  ASSUMED  TO  BE  PURE  REAL 

IMPLICIT  COMPLEX  (C)  j 

REAL  IHYP* IHYPO  l 

ARG (CZZ ) *ATAN2 ( A1MAG (CZZ)iREAL(CZZ) ) 

SET  UPPER  LIMIT  ON  NUMBER  OF  TERMS 
LMT =500 
ARZG=ARG(-CZ) 

ARGZ=ARG(CZ) 

1000  FORMAT ( 'ONC  CONVERGENCE* / *  CHYP* * , 2( IPE10.3,2X) , 5X, *CHYPO* • ,2 ( ElOj 

13, 2X))  1 

1001  FORMAT ( *  I*  * , I5,5X, 'ALPHA*  1 , F10. 5, 5X, ' GAMMA*  • .F10.5.5X, *CZ*  ',2 
1(1PE10.3,2X),3X,*RHYP>  * , ElO. 3,5X, • IHYP*  *,E10.3> 

300  CHYP0=<1.,0. )  ! 

1*1 

Z»REAL(CZ) 

TEST  FOR  LARGE  ARGUMENT 

IF(CABS(CZ).GT.3.)G0  TO  500 

SMALL  ARGUMENT  EXPANSION 

IHYP*0. 

RHYPO*ALPHA/GAMMA*Z 

SRHYP=RHYPO  | 

320  1*1*1 

CHECK  NUMBER  OF  TERMS  t 

IF (  I.GT.LMT )G0  TO  340 
K  =  I-1 

RHYP*RHYPO* ( ALPHA+K ) / (GAMMA+K )*Z/FLOAT ( I ) 

SRHYP=SRHYP*RHYP 
TEST  FOR  CONVERGENCE 

I F (  ABS ( RHYP— RHYPO )  «LT • l.E-05 )G0  TO  330 
RHYPO=RHYP 

GO  TO  320  1, 

WRITE  ERROR  MESSAGE 
340  WR I TE ( 6 , 1000IRHYP, RHYPO, CHYPO 

WRITE(6,1001) I, ALPHA, GAMMA, CZfRHYP, IHYP 
330  RHYP=SRHYP*1. 

RETURN 

500  CONTINUE 

ASYMPTOTIC  EXPANSION  FOR  LARGE  ARGUMENT 


RSUM~1. 
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S SUM* la 
RTERMO-1. 

STERMO-1. 

CALL  GAMMB ( ALPHA* GAMMA# A* B* BMA ) 

C  SET  UPPER  LIMIT  ON  NUMBER  OF  TERMS 
LMT-5 

00  600  I-ltLMT 
K*I-l 

RTERM*<ALPHA+K)*(ALPHA-GAMMA+l*K)/tFLOAT( !>•(  -2 ) )*RTERM0 

RSUM-RSUM+RTERM 

RTERMO-RTERM 

STERM* ( GAMMA- ALPHA+K ) * ( 1 .-ALPHA-** ) / < FLOAT ( I ) *  Z ) * STERMO 
SSUM*SSUM«-STERM 
STERMO- STERM 
600  CONTINUE 

CF*CABS IC2 ) ** ( -ALPHA ) *CEXP  U  0. • 1 . )* ( -ALPHA ) *ARZG ) /BMA*RSUM 
IF(CABS(CZ).GT.99. )G0  TO  700 

C  F-CF +C ABS ( CZ ) ** I ALPHA-GAMMA ) *CEXP ( 1 0. * 1 . ) * I ALPHA-GAMMA ) *ARGZ +CZ ) / 
IA*SSUM 
700  CONTINUE 

RHYP*REAHCF)*B 
I  HYP* A I MAG ( CF ) *8 
RETURN 
ENO 


SUBROUTINE  GAMMB(AVB« A1»B1.BMA) 

THIS  SUBROUTINE  CALCULATES  GAMMA  FUNCTIONS  FOR  CONFLUENT  HYPERGE0METR1C 
FUNCTION'S  LARGE  EXPANSIONS 

B1*GAMMA(B) 

I F ( A.GT.B) GO  TO  200 
BMA*GAMMA I B-A ) 

ISO  CONTINUE 

IFtA.LT.O.IGO  TO  100 
A1*GAMMA( A) 

RETURN 

100  Al*-3. 141592/ (GAMMA < -A+l . )*S INI-3. 14159*A)J 

RETURN 

200  BMA*-3. 161592/ (GAMMA  ( A-B+l  •  »*SIN(  3. 141592*U-8m 

GO  TO  150 
ENO 
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Col umns 


1-10 

11-20 

21-30 


RHOS 

DELRH 

XNR 


the  first  value  of  p 

the  increment  of  p 

the  total  number  of  increments. 


Again  the  variables  are  punched  in  ten  column  fields  in  decimal  form. 

If  XNR  =  0,  then  the  program  reads  a  new  parameter  card. 

(b).  The  Von  Karmann  spectrum  calculations  of  the  mean  square 
angle  of  arrival,  two  point  correlation  functions,  and  structure 
functions  mentioned  earlier  are  performed  by  numerical  integration.  The 
programs  for  both  plane  and  spherical  wave  input  are  shown  in  listings 
(B-2)  and  (B-3) ,  respectively.  The  integrations  are  performed  by 
using  a  96  point  Gaussian  quadrature  algorithm.  The  Bessel  functions 
were  evaluated  using  a  polynomial  expansion  from  Handbook  of  Mathematical 
Functions  by  Abramowitz  and  Stegun. 

The  weighting  factors  and  arguments  for  the  Gaussian  inte¬ 
gration  are  read  as  input  data  and  are  shown  in  listing  (B-4).  All 
user  input  data  should  follow  these  cards.  Since  both  programs  use 
the  same  input  data  only  one  data  deck  will  be  explained. 

The  first  user  input  data  card  specifies  the  parameters  for 
the  angle  of  arrival  correlation  function  calculation.  The  data  is 
again  punched  in  ten  column  fields  in  decimal  form.  The  input  variables 


are: 

Columns 

1-10 

XL0  = 

t0,  the  inner  scale  in  meters 

11-20 

XL  = 

L,  the  range  in  meters 

21-30 

XLAMB  = 

x,  the  wavelength  in  meters 

31-40 

XCAPL0  = 

L0,  the  outer  scale  in  meters 

41-50 

A  = 

A,  the  outer  scale  normalizing  constant. 

If  XL  =  0,  the  correlation  function  calculations  are  terminated  and  the 
program  begins  the  structure  function  calculations. 

Otherwise,  the  program  reads  the  separation  parameters:  starting 
value,  increment,  and  total  number  of  increments.  These  variables  are 
punched  in  ten  column  fields  in  decimal  form.  They  are 
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LISTING  B-2 


C  ****  IC.6  M ICKCN  PRUPAGA f ION  PHENOMENA  PROGRAM  ♦*** 

C  *****  ANGLE  OF  ARRIVAL  AND  RELATED  STRUCTURE  FUNCTION  CALCULATION 
C  FLR  THE  VON  KARMANN  SPECTRUM*  ANO  PLANE  INPUT  **♦*♦ 

EXTERNAL  FC T « FQ! 0* FC 1  I • FCT2, FCT 3 

KCAL*fa  YINT«YINTG,FCT,YL,<U*YINTl,YINT2*XKAP 

COMMON  RHQ*XL«XK«XKMfXKQ 

COMMON  /INTG/TEST 

CALL  GAUSZ 

TFS  T =0 • 

C  READ  INPUT  DATA 

IOC  REAU(5*900)XLC,XL«XLAMB,XCAPL0*A 
SCO  FORMAT (SEIC.3) 

IFIXL.LC.O. )GC  TO  20i 
C  CALCULATE  CONSTANTS 
XKM*S. 94/XLO 
XKL*A/XCAPLO 
PI  =  3. 1415927 
XK* 6. 2831 84/XL AMB 
C  READ  ARGUMENT  VALUES 
ISO  READ ( S * 900 ) RHOS*  DELRHO»XNR 
IF  (XNR.EC.O. )G0  TO  100 
C  CALCULATE  ANGLE  OF  ARRIVAL  CORRELATION 
NR-XNR 

NR  I IE ( 6* 1002) 

1002  FORMAT  1 1 1 ANGLE  OF  ARRIVAL  CORRELATION  FUNCTION  •) 

WRITE(6«100S ) XL*  XLO*  XCAPLO*  A*  XLAMB 

100S  FORMAT!1  RANGE* • * IPE 11 . 3, •  M' »SX* • INNER  SCALE*  ',E10.3*»  M»,5X, 

1 1 OU TER  SCALE*  *«E10.3**  M»,5X,»A*  » , E 10. 3, 5X * ' LAMBDA*  SE10.3*'  M* 
2) 

DO  200  1*1 »NR 

RHC=RHGS*DELRHC*FLCAT( I— l I 
C  ELEVATION  ANGLE  OF  ARRIVAL  CORRELATION 
CALL  I NTGRL (0.*l.E-3*8«FCT*YINT) 

BALPHA=2.*PI**2*XL*YINT 
IF  I RHG.EC.O. ) B AN* B ALPHA 

balphn=balpha/ban 

WRITE  1 6, 10G0 ) RHQ * BAL PHA *  BALPHN 

1000  FORMAT! *ORHO=  '* IPE 10.3* SX»  *B ALPHA*  '  , E10 . 3, 5X, ' BALPHAN*  *,E10.3) 

C  AZIMUTH  ANGLE  OF  ARRIVAL  CORRELATION 

CALL  INTGRL(0.,l.fc-3,8,FCT0,YINT) 

BBETA=2.API**2*XL*YINT 
I F ( KHG.EO.O. ) B BN* OB ETA 
B8ETAN*BBETA/B8N 

1001  F0RMAT(21X* 'BBETA*  » , 1  PE  10. 3 ,SX , • BBETAN*  ',E10.3) 

200  WRITE(6«1001)6BETA*BBETAN 

GO  TO  ISO 

C  CALCULATE  STRUCTURE  FUNCTIONS 
C  READ  INPUT  DATA 
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20 1  REALMS,  9C0>XL(J,XL»XLAMB,XCAPLG,A 
IFIXL.EC.O. )CALL  EXIT 
C  CALCULATE  CONST  ANTS 

XK=6.2b3l484/XLAM8 
XKM  =  5.<)4/XL0 
XKC=A/XCAPLU 
WKITHbt  1003) 

1003  FORMAT ( * ^STRUCTURE  FUNCTIONS  AND  SPECTRA  •) 
WKI TE ( 6* 1005 ) XL , XLO, XCAPLO, A , XLAMB 
C  Rt  AO  ARGUMENT  VALUES 
230  REA0(5f900)RtlGStDELRHCfXNR 
IFIXNR.EC.O. )GC  TO  201 
NR=XNR+ l • 

00  220  1=1 « NR 

RHC=RHCS  +  OFLRHO*FLC A T (  I-l) 


XKAP=RHC  p 

C  REFRACTIVE  INDEX  FLUCTUATION  SPECTRUM 
PHlK=PHlFPS<XKAP,XKM,XKO> 

C  REFRACTIVE  INDEX  FLUCTUATION  STRUCTURE  FUNCTION 
CALL  INTGRL(1«E”03»1«E-02»8»FCT1»YINT1) 

0N=8.*PI*YINT1  *» 


C  PhASL  STRUCTURE  FUNCTION 

CALL  INTGRU  1  .  E-02, 1 .0E-01 ,6,  FCT2 1 Y  INT2  ) 
US=4«*PI**2*XK**2*XL*Y  INT2 


C  WAVE  STRUCTURE  FUNCTION 

CALL  INTGRLI l . E-03, 1 . E-02, 7 , FCT3, Y INT2 ) 

DSW=8.*PI**?*XK**2*XL*YINT2 
220  WRI  TE(6, L004)RHO,PHIK,0N,OS*DSW 

1004  FORMAT ( *OARGUMENT=  • , 1PE 10 . 3 * / , •  REFRACTIVE  INDEX  FLUCTUATUION  SPE 
1C  TROM*  *,E10. 3, 5X, 'STRUCTURE  FUNCTION*  ',Ei0.3,/,»  PHASE  STRUCTURE 
2  FUNCTION*  ,»E10.3,5X*' WAVE  STRUCTURE  FUNCTION*  SE10.3) 

GO  TO  230 
END 


0 

0 

0 


c 


DOODLE  PRECISION  FUNCTION  FCT(XKAP) 

INTEGRAND  FOR  CALCULATING  ELEVATION  ANGLE  OF  ARRIVAL  CORRELATION 
REAL+8  ARGl ,FCT, XKAP 
COMMON  RHOtXLtXKyXKMfXKO 
ARG3=(XKAP/XKM)**2 
IF (AKG3.GT.174. )GOTO  ICO 
ARG2*XK AP*RHO 
ARGl=XKAP**2*XL/XK 
IF ( ARGl. EC. 0. >GU  IC  100 

FCT=( l.+OS INI ARGl ) /ARG 1) * ( BESO ( ARG2 )-BES IX( ARG2 ) )*XKAP**3*PHIEPS 
IKAP.XKM,XK0) 

RETURN 


0 

0 

0 


100  FCT=0. 
RETURN 
END 


y 


i  J 


L 
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DUUULF  PRFC I  SIGN  FUNCTION  FCTO(XKAP) 

C  INIFGRAND  FCK  CALCULATING  AZIMUTH  ANGLE  OF  ARRIVAL 
REAL*8  AKGlfFCTO* XKAP 
COMMON  RHt1»XLfXKfXKMyXKO 
ARG3= ( XKAP/XKM ) **2 
IFIARG3.GT.174. )GCTU  100 
ARG2-XKAP*RH0 
ARGl=XKAP**2*XL/XK 
IF (ARGl.FC.O. IGO  TO  100 

FCTC-I  l.+OSINI  ARGl )/ ARGl )*BFS IX (ARG2)4XKAP**34PHIEPS( XKAP, XKM.XKO) 
RETURN 
100  FC  T0=0. 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  FCTKXKAP) 

C  INI EGRANL)  FOR  CALCULATING  INDEX  STRUCTURE  FUNCTION 
REAL*b  ARG? *FCTl»XKAP 
COMMON  RHCyXLyXK , XKM, XKQ 
ARGl*XK AP**2*XL/XK 
ARG2=XKAP*RH0 
ARC  3*1 XKAP/XKM 1 4*2 
I F ( ARG3.G  T .174.1GOTO  100 
1FIARGI.E0.0.IGU  TO  100 

FCI l*( 1.-USIN(ARG2)/ARG2)*PHIEPS(XKAP,XKM,XK0)*XKAP**2 
RETURN 
100  FCTI*0. 

RETURN 

END 


OOUBLE  PRECISION  FUNCTION  FCT2IXKAP I 
C  INI EGRAND  FOR  CALCULATING  PHASE  STRUCTURE  FUNCTION 
REAL*8  ARG1*FCT2»XKAP 
COMMON  RHOyXLvXKyXKMyXKO 
ARG 3= (XKAP/XKM)** 2 
IFIARG3.GT.174.) GO  10  100 
ARG1*XKAP**2»XL/XK 
ARG2=XKAP*RH0 
IFIARGl.EQ.O. )G0  TO  100 

FCT2=(l.-BES0( ARG2) ) * ( 1 . *DS IN ( ARGl ) /ARGl ) *PHI EPS C XKAP f XKM.XKO ) *XKA 
IP 

RETURN 
100  FCT2=0. 

RETURN 

END 
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UCLBLE  PRECISION  FUNCTION  FCT3(XKAP) 

C  1NICGRANO  FOR  CALCULATING  THE  WAVE  STRUCTURE  FUNCTION 
KF.AL*H  XK AP  t  FCT  3 
CCRFUN  RHG'XLtXK'XKN'XKO 
ARG2=XKAP*RHG 
ARG3=(XKAP/XKM)**2 
IF<ARG3.GT.174.)GUTO  LOO 

FCI 3=  ( 1.-BESGIARG2) )*PHIEPS(XKAP*XKM* XKU)*XKAP 
RETURN 
ICO  FCT  3=0. 

RE  I  URN 
END 


FUNCI1CN  PFtltPSIXK  »XKM»XKO  ) 

C  REFRACTIVE  INUEX  FLUCTUATION  SPECTRUM 
RFAL*8  XK 
ltSTl=(XK/XKM)**2 
l F ( TESTi.GT.l 7A.  )GC  TO  LOO 
A  =  EXP ( -TtSTi ) 

b=(XK**2*XK0**2)**(ll./6. ) 
IFST2=ALCG10(A)-ALGG 10(B) 

IF<  ABSM  EST2I.GT.60.  )GC  TO  LOO 

PHItPS=.C33*A/e 

RETURN 

LOO  PH  I EPS  =  C. 

RE  I URN 
END 
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SUbRlMJ!  i (ML  INTGRLIXLS, XUS, NINT,FCT, VINT) 

C  URuEK  OF  MAGNITUDE  INTLGRAT ION  ROUT  INF 
COMMON  /1MTG/TEST 
REALMS  YlNT,YlNTG, FCT,XL,XU 
EXTERNAL  FCT 

1G03  FURMATI'  I  NT£GRAL=  *  »  IPE 10. 3  ♦  *>X ,  •  LOWER  LIMIT=  »,ei0.3,*>X, 
1  *  UPPER  LIMIT*  »,EI0.3) 

100?  FORMAT!//) 

IF! TEST. GT.O. ) WRITE  1 6,1002) 

XL=XLS 

XU=XUS 

CALL  DUG96 1  XL ♦ XU  »FCT»YlNT) 

IFITEST.GT.O. )WR1TEI6, 1003 ) Y I  NT, XL  *  XU 

UU  300  J=t,NINT 

XL=XU 

XU=XU*10. 

CALL  UCG96 I  XL  »  XU » FCT , Y INTG) 

YINT=YINT+YINTG 

IFITEST.GT.O. ) WRITE  1 6, 1003 )YINTG, XL, XU 
300  CONTINUE 
RETURN 
ENO 


SUUKGU TINE  DGG96(XL»XU»FCT »YINf ) 

C  96  POINT  GAUSSIAN  QUADRATURE  INTEGRATION  ROUTINE 
IMPLICIT  RFAL*8  (A-H,0-Z) 

COMMON  /GAUSS/X(48),W(48) 

SUM=( XU+XL ) /2. 

U  IF* I XU-XL ) /2. 

YiNT*0. 

UC  100  1*1,48 
ARGl=OIF*XII) ♦SUM 
ARG2=”U IF*X 1 1 ) ♦SUM 

ICO  YINT=Wt I ) * ( FCT (ARG1 1+FCTIARG?)  )+YINT 
Y INT=YINT*0IF 
RETURN 
ENC 
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SUBRUU  I  l  M  GAUSZ  I 

C  KlAU  AM  SlOKfc  GAUSSIAN  INT  EGKAT  IGN  ARGUMENTS  AND  WEIGHTS  I 

IMPLICIT  HE  AL*d  ( A-H»C-Z  I 

CGMRON  /GAUSS/X(4ti)«W(4b)  j 

KEAL(6,999)X 

HEAL (6,999 )W 

999  F0RMAH5E16.15)  r 

MRITEIbtllillX 
MR1TE(6«111I)W 

1UI  FORMAT  (IX«6(F17.15,3X)/) 

RETURN 

LNC  I 


\] 

FUNCTION  BESG(X) 

KE AL*8  THETA  r 

C  JOIXI  I 

IFlX.GT.a. )GC  TO  IOO  *' 

Y=X/3. 

BESC=I  .0-2.24999*Y**2U.26562*Y**4-.3l638  7*Y**6+.044447*Y**8-3*944' 
14E-3*Y**10*2.10000E-4*Y**12  l 

RETURN 

10U  Y=3 ,/X  r 

ThETA=X-.78639B-.04l663*Y-3.9 >4E-05*Y**2*2.62573E-3*Y**3-5.4125E-4 
l*Y**4-?.9333E-4*Y**5+l.3558E-4*Y**6 
FG=.79/8B46-7.7E-7*Y-5.5274E-3*Y**2-9.512E-5*Y**3+1.37237E-3*Y**4- 
17.2805E-4*Y**5+1.4476E-4*Y**6 

BESC=EC*DCGS(THETA)/ SORT  C  X )  * 

RETURN 

ENL  f] 


FUNG  I  I CN  BEST  X  (  X  I  J_j 

KEAL*8  THETA 

C  JUXJ/X  r] 

IFIX.GT. 3. )GG  TO  IOO 
Y=X/3. 

bESlX=.50-.56?499*Y**2*.2l0935*Y**4-3.95428E-2*Y**6*4.433l9E-3*Y** 
18-3.1761C-4*Y**10+1.109E-5*Y**12  jl 

RETURN  L 

100  Y=3./X 

THE TA=X-2. 356 19*. 124996* Y*5.650E-05*Y**2-6,37879E-3*Y**3+7,4348E  41 
l*Y**4+/.9824fc-4*Y**5-2.9l66E-4*Y**6 

F 1=.797884+1.56E-6*Y+1.65966E-2*Y**2-1.7105E-4*Y**3-2.49511E-3*Y**' 
14+1.13653E-3*Y**5-2.0033E-4*Y**6 
bESlX  =  Fl*OCOS( THETA) /(X*SQRT(X)  ) 

RETURN  ^ 

END 


LISTING  B-3 


C  ****  10.6  M ICKCN  PRUPAGAT  ION  PHFNOMENA  PROGRAM  **** 

C  *****  ANGLE  OF  ARRIVAL  ANU  RELATED  STRUCTURE  FUNCTION  CALCULATION 
C  FOR  THE  VON  KARMANN  SPCCTRUM,  AND  SPHERICAL  WAVE  INPUT  ***** 

EXTERNAL  FCT  *  F  Cl  ti ,  FCT 1 »  PC  I  2 , FC  T  3 
REAL*b  YINT,YINTG,ECT,YL,XU»Y1NTI,YINT2,XKAP 
COMMON  RHUtXLt XK.XKM.XKO 
COMMON  /INTG/TEST 
CALL  GAUSZ 
TfcST=0. 

C  KF AD  INPUT  DATA 

100  READ(5,900)XLO,XL,XLAMB,XCAPLO,A 
900  FORMAT (5E10. 3) 

IF (XL.EG.O. )G0  TO  201 
C  CALCULATE  CONSTANTS 
XKM=5. 94/XLG 
XKC*A/XCAPLO 
PI *3.1  A 159? 7 
XK=6.28318A/XLAMB 
C  READ  ARGUMENT  VALUES 
150  REAU(5,900)RHUSf DFLRHO,XNR 
I F ( XNR • EC • 0. ) GC  TO  100 
C  CALCULATE  ANGLE  OF  ARRIVAL  CORRELATION 
XNR=XNR-l. 

NR=XNR 

WR  I  TE ( 6 , 1002  1 

1002  FORMAT  ( MANGLE  OF  ARRIVAL  CORRELATION  FUNCTION  •) 

WRITE (6tl005)XL*XLUt XCAPLQ, A, XLAMB 

1 005  FORMAT!1  RANGE*'  ,  IPL11.3,  •  M',5X, 'INNER  SCALE*  ',E10.3,'  M',5X, 

I'OUTER  SCALE*  '.E10.3,'  M',5X,'A=  ' ,L 10. 3,5X, 'LAMBDA*  SCIQ.3,'  M' 
?) 

DO  200  1*1, NR 

KHC=RHCS*UELRHG*FLCAT (1-1) 

C  ELEVATICN  ANGLE  OF  ARRIVAL  CORRELATION 
CALL  IN!GRL(0.,l.E-3,8,FCr,YINT) 

BALPHA=.5*PI**2*8.*Y I  NT 

IF ( RHO. EC. 0. ) BAN* 8 ALPHA 

BALPHN=B ALPHA/ BAN 

WRITE (6, IOOC)KHO,BALPHA,BALPHN 

1000  FORMAT! «ORHO=  ' , 1PE10. 3,5X, * BALPHA*  ' ,E10.3,5X, 'BALPHAN*  ',E10.3) 

C  AZIMUTH  ANGLE  OF  ARRIVAL  CORRELATION 

CALL  INTGRL(0.,1.F-3,8,FCT0,Y1NT) 

BBETA=.5*Pl**2*8.*YINT 
IF(RH0.EC.0.)8BN=BBETA 
b BE  T  AN* 8 BET A/ CBN 

1001  FORMAT  (21X,(’BBETA*  '  ,  l PE  10.  3 , 5X,  ' 88ETAN*  SE10.3) 

200  WR 1TE(6, 1001)B BE TA,B6ET AN 

GC  TO  150 

C  CALCULATE  STRUCTURE  FUNCTIONS 
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C  READ  INPUT  DATA 

201  REAU(t>»900)  XLO*XL  *XL AMB*  XCAPLO*  A 

1F(XL.EQ.0.)CALL  EXIT 
C  CALCULATE  CONSTANTS 

XK=6.2B31484/XLAMB 
XKM=S. 94/XLO 
XKO=A/XCAPLO 
WRITEI6.1003) 

1003  FURMAT ( • ISTRUCTURE  FUNCTIONS  AND  SPECTRA  • ) 

C  READ  ArtGUMENI  VALUES 

230  READ!  .>*900) RHUS*  DELRHO*  XNR 
IFIXNR.ECl.O.JGC  TO  2 01 
NR=XNR+1. 

00  220  1*1 « NR 
RH0=KH0S+DELRH0*FL0ATII-1) 

XKAP=RHO 

C  REFRACTIVE  INDEX  FLUCTUATION  SPECTRUM 
PHIK=PHlEPS(XKAP»XKM,XKO> 

C  REFRACTIVE  INDEX  FLUCTUATION  STRUCTURE  FUNCTION 
CALL  INTGRL(l.E-03,l.E-02. 8, FCTl » YINT 1 ) 

DN=8.*PI*YINTl 
C  PHASE  STRUCTURE  FUNCTION 

CALL  INTGRL(l.E-02,1.0E-0l*6,FCT2,YlNT2) 

DS=YINT2*PI**2*XK**2*8. 

220  WRITE (6*1004 ) RHO*  PH  I K  *  DN*  DS 

1004  FORMAT ( ' 0 ARGUMENT*  • , IPEIO. 3,/, •  REFRACTIVE  INDEX  FLUCTUATION  SP 
ICTRUM*  *,E10. 3, 5X, ‘STRUCTURE  FUNCTION*  ‘.ElO.*,/**  PHASE  STRUCTUR 
2FUNCT 1  ON*  ' ,E 10. 3 ) 

GO  TO  230 
END 


OCCBLE  PRECISION  FUNCIION  FCT(XKAPPA) 

C  IN IEGRANO  FOR  CALCULATING  ELEVATION  ANGLE  OF  ARRIVAL  CORRELATION 
EXTERNAL  FCT4 

RFAL*8  ARG1.FCT , XK AP , CXL * DSPHR , XK APP A , FCT4 

COMMON  RHO,XL,XK,XKM,XKO 

COMMON  /SPHERE /XKAP 

XKAP=XK APPA 

UXL-XL 

ARG2=XKAP*RH0 

CALL  DCGb  IO.,OXL,FCT4*DSPHR) 

FCT=IBESC< ARG2)-BES1X(ARG2) )*XKAP**3*DSPHR 

RETURN 

END 
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UCUI3LE  PRECISION  FUNCTION  FCTO(XKAPPA) 

C  INTEGRAND  FOR  CALCULATING  AZIMUTH  ANGLE  OF  ARRIVAL 
EXTERNAL  FCT4 

RE A L*B  ARG1 , FC TO , XK AP , OXL » OSPHR , XK APP A , FC T4 

COMMON  RF(],XLtXK,XKM,XKO 

COMMON  /SPHERF/XKAP 

XKAP=XK APP A 

ARG2=XKAP*RHC 

DXL-XL 

GALL  DCG«  (0. ,DXL,FCT4,DSPHR) 

FU C=bFS  IX  (  AKG2)  *XKAP** 3*USPHR 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  FCTUXKAP) 

C  INTEGRAND  FOR  CALCULATING  INDEX  STRUCTURE  FUNCTION 
REAL*B  ARG2 *FC T I »  XKAP 
COMMON  RHO,XLtXK,XKM,XKU 
ARGl=XKAP**2*XL/XK 
ARG2=XKAP*RH0 
ARG3=(XKAP/XKM)**2 
IF( ARG3.GI . 17'». IGGTO  100 
IF (ARG1.LG.0. )G0  10  LOO 

FCT I = ( l.-DSINl ARG2)/ARG2)*PHIEPS(XKAP,XKMtXK0)*XKAP**2 
RETURN 
ICO  FC  T I  =  0. 

RETURN 

LND 


DOUBLE  PRECISION  FUNCTION  FCT 2 ( XKAPPA ) 

C  INTEGRAND  FOR  CALCULATING  PHASE  STRUCTURE  FUNCTION 
EXTERNAL  FCT4 

REAL*M  ARGltFCT?,XKAP,DSPHR,UXL,XKAPPA,FCT4 

COMMON  RHO,XLtXK,XKM,XKO 

COMMON  /SPHERt/XKAP 

XKAP=XK APPA 

ARG2=XKAP*RH0 

UXL=XL 

CALL  UCGH  (0.,CXL,FCT4,DSPHR) 

FCT2=( l.-OLSOI ARG2) )*DSPHR*XKAP 

RETURN 

END 
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uCLULE  PRCCISICN  FUNCTION  FCT3TXKAP) 

C  INTEGRAND  EUR  CALCUL A T ING  THE  WAVE  STRUCTURE  FUNCTION 
KfcAL*b  XKAP*  FC  T  3 
COMMON  RhOt XLtXK.XKMt XKU 
ARG2=XKAP*RHG 
ARG3=(XKAP/XKM)**2 
IN  AKG3.GT.  17*.  )GGTU  100 

FCTi=( l.-BESGT ARG2) )*PHIFPS (XKAP,XKMtXKO>*XKAP 
RETURN 
ICO  FCT  3  =  0. 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  FCT  4TDETA ) 

REAL+8  XKAP*DETA*FCTA 
COMMON  RHO  tXLtXKfXKMtXKO 
COMMON  /SPHERF/XKAP 
ARGl=XKAP*Xl/OETA 

FCTA-=(XL/0ETA+CC0S(XL*(XL-DETA)  +  XRAP**2/(  2.*XK*DETA )  )  )  **2*PHI  EPS  I 
1 ARG 1 »  XKM  »  XKO ) 

RETURN 

ENC 


FUNCIICN  PHIFPS(XK»XKM» XKO ) 

C  REFRACTIVE  INOEX  FLUCTUATION  SPECTRUM 
REAL*fa  XK 
TESTl=(XK/XKM)**2 
iriTESTl.GT.17A. )GG  TO  100 
A*EXPI-TEST1) 

B=(XK**2*XK0**2)**(ll./6. ) 
TEST2-ALCG10TAJ— ALCG10TBT 
IF(ABS(TEST2).GT.60. ) GO  TO  100 
PHIEPS*.033*A/E 
RETURN 

ICO  PHIEPS=0. 

RETURN 

ENC 
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SUbKOUf  INE  INTGRLlXLS, XUS , N l NT , FC T , Y INT ) 

C  URUER  CF  MAGNlTUUF  l NTEGRAT ION  ROUT INE 
COMMON  /INTG/1FST 
REAL*8  YlNT,YlNTG,FCT,XL,XU 
EXTCRNAL  FCT 

LOG3  FORMAT ( '  INTEGRAL*  '  ,  IPE 10 . 3, 5X , • LOWER  LIMIT*  »,El0.3,5X, 
1  *  UPPER  LIMIT*  »,E10.3) 

1002  FORMAT!//) 

IF ( TEST. GT. 0. ) WRITE (6, 1002) 

XL=XLS 

XU=XUS 

CALL  0CG06 ( XL , XU, FCT » Y I  NT ) 

IF! TEST.GT.O. ) WR ITE ! 6 , 1003 ) Y I  NT , XL ♦ XU 
00  300  J= It  MINT 
XL  =  XU 
XU=XU*10. 

CALL  DCG96!XL,XU,FCT,YINTG> 

Y INT*Y INT+YINTG 

IF(TEST.oT.O. ) WRITE! 6, 1003 ) Y INTG, XL, XU 
300  CONTINUE 
RETURN 
ENU 


SUBROUTINE  DCG96 { XL , XU , FCT , Y I N T  ) 

C  96  PUINT  GAUSSIAN  CUADRATURF  INTEGRATION  ROUTINE 
IMPLICIT  KEAL*8  !A-H,0-Z) 

COMMON  /GAU$S/X!48),W!48) 

SUM*! XU+XL ) /2. 

U IF*! XU-XL  )  /2 • 

YINT=0. 

00  100  1*1,48 
ARGl*OIF*X!  D+SUM 
ARG2=-UIF*X!  I  HSUM 

100  Y INT*W ( I ) * ( FCT ( ARGl ) +i  C •  > ARG2  )  )  ♦Y1NT 
Y INT=YINT*CIF 
RETURN 
ENU 


SUBROUTINE  GAUSZ 

C  RE AO  ANO  STCRE  GAUSSIAN  INTEGRATION  ARGUMENTS  AND  WEIGHTS 
IMPLICIT  REAL*8  ( A-HtC-Z ) 

COMMON  /GAUSS/X!48),W!4H) 

REA0!S,S99)X 
REAC !  b ,999  )  W 
999  FORMAT ! 5F16. lb ) 

WRITEC6,  1 1 1 1 )  X 
WRlTE(6,llll)W 

llll  FORMAT (IX,6(F17.  IS, 3X)/) 

RETURN 

END 
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[ 

FUNGI  1 CN  BE  SC ( X ) 

C  JC  <  X I 

1 F ( X.C  T • 3. )Gl!  TO  100 
Y=X/3. 

8ESI.M  .0-2.24999*Y**2*l.26S62*Y**4-.3l638  7*Y**6+.044447*Y**8-3.944 
l4E-3*Y*+lO+2. 10000fc-4*Y**l2 

RETURN  1 

100  Y=3./X 

lHETA  =  X-.78b398-.04L663*  Y-3.9S4E-0*)*Y**2  +  2.62573E-3*Y**3-S.4125E-4 
1*Y**4-2.93  33E-4*Y**S  +  I.3t»b8fc-4*Y**6  I 

FC  =  .7978b46-7.7C-/*Y-5.52  74E-3*Y**2-9.‘»126-5*Y**3^1.3723/E-3*Y**4- 
!7.2805E-4*Y**b+1 .4476E-4*Y**6 
BbSC=FC*CGStlHFTA)/SC)RT (X) 

RE  1  URN  * 

fcNC 


FUNCTION  BES1X ( X ) 

KEAL+H  THETA 

C  JKXJ/X  f 

IF(X.GT.3. )G0  TO  100 

Y=X/3.  1 

HbSiX=.50-.bo2499»Y**2*.2L0935*Y**4-3.95428E-2*Y**6+4.433l9E-3*Y** 
lb- 3.  I 76lb-4+Y**10+l. 109E-S*Y**12  f 

RETURN  I 

100  Y  =  3 • / X 

THETA=X-2.356t9*.l24996*Y*l>.6b0E-05*Y**2-6.37879E-3*Y**3*7.4348E-4j 
1*Y**4  +  7.9824E-4*Y**5-2.9166E-4*Y’**6 

Fl  =  ./978b4M.56E-6*YM.65966F-2*Y**2-l.71C5E-4*Y**3-2.495llE-3*Y** 

14+1.13653E-3*Y**5-2.0033E-4*Y**6 

BES1X=F1+0C0SUHETA)/IX*SGRT(X)  )  j 

RETURN  I 

END 

l 

i 

1 


' 
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Col umns 


1-10  RHOS  =  starting  value  of  p 
11-20  DELRHO  =  increment 

21-30  XNR  =  total  number  of  increments. 


If  XNR  =  0,  a  new  parameter  card  is  read,  otherwise  another  set  of 
separation  parameters  is  expected. 

The  data  deck  for  the  structure  function  calculations  is 
identical  to  the  correlation  data  deck.  The  same  variables  are  used 
and  the  same  action  is  taken  for  XNR  =  0.  When  XL  =  0,  the  program 
is  terminated. 


2.  Large  Aperture  Calculations 

Two  programs  were  written  to  perform  the  calculation  of  the 
large  aperture  results.  One  calculates  the  large  aperture  angle  of 
arrival  and  linear  phase  error;  the  other  calculates  the  error  between 
the  large  and  small  aperture  results.  Both  programs  use  the  Von  Karmann 
spectrum. 

The  large  aperture  angle  of  arrival  program  shown  in  listing 
(B5)  performs  four  nested  integrations  to  calculate  <ox,  <a£>  and 
<a2k>.  Gaussian  quadrature  integration  of  different  orders  is  used  in 
each  case.  The  Gaussian  subroutines  are  contained  in  the  IBM  Scientific 
Subroutine  Package.  They  are  essentially  the  same  as  DQG96  used  earlier 
except  the  weighting  factors  and  arguments,  are  contained  in  the 
subprograms.  The  large  aperture  error  program  evaluates  the  error 
function  a(<R)  in  a  similar  manner. 

The  input  data  for  both  programs  is  read  with  a  namelist  format. 
The  variables  are 


XCAPL0  =  L0,  the  outer  scale  in  meters 
XL0  =  Hq  ,  the  inner  scale  in  meters 
RCAP  =  R,  the  aperture  radius  in  meters 

A  =  A,  normalizing  constant  in  Von  Karmann  spectrum 
XLAMB  =  x,  wavelength  in  meters 
XL  =  L,  range  in  meters. 


After  each  calculation  a  new  data  card  is  read.  The  program  is 
terminated  by  setting  XLAMB  =0. 


LISTING  B-5 


C  ****  10.6  MICRON  PROPAGATION  PHENOMENA  PROGRAM  •*** 

C  ****  ANGLE  OF  ARRIVAL  AND  ERROR  CALCULATION 
C  USING  THE  PHASE  STRUCTURE  FUNCTION  **♦* 

IMPLICIT  REAL*8  (D) 

COMMON  R,RHO« GAMMA, RCAP 

COMMON  /PHI/XKMtXKO 

COMMON  /PHASE/XK,XL 

COMMON  /APRAD/DOCAP 

EXTERNAL  OINTl,OOEN, OERFCT 

NAMEL I  ST/ INPUT/XCAPLOf  XLO, RCAP , A, XLAMB.XL 

PI»3. 1415927 

WRI TE (6*1000) 

1000  FORMAT ( ' I ANGLE  OF  ARRIVAL  CALCULATION  USING  PHASE  STRUCTURE  FUNCTI 
ION*/*  (RESULTS  NORMALIZED  TO  CN**2 )  •  ) 

1005  FORMAT ( '  RANGE3*  ,1PE11.3«*  M*  »5X»  *  INNER  SCALE3  *,El0.3t*  M',5X, 

I  *  OUTER  SCALE3  *,E10.3,'  M*,5X,*A*  *, ElO. 3, 5X, 'LAMBDA3  SE10.3,*  M* 

2) 

100  READ(5f INPUT) 

WRI TE(6t 1005) XLtXLOtXCAPLOf A,XLAMB 

I F ( XLAMB.EQ.O. )CALL  EXIT 

XKS6.28318/XLAMB 

XKM»5. 92/XLO 

DRCAP*RCAP 

XKO=A/XCAPLO 

CALLDQGI2(0.D0f DRCAPt DINTltDRINTl ) 

CALL00GI2 (O.DOf DRCAP tDDENtORDEN) 

GX0*-4.*PI ♦*2*XX**2*XL*DRINTl/ ( PI*DRDEN**2 ) 

A0A*GX0/XK**2 

DDCAP*2.*DRCAP 

CALL  0QGI2(0.f DDCAPf DERFCTtDRNTl ) 

DERR0R34.*PI*XK**2*XL*0RNTI/DRCAP**2 

DPSERR3DERR0R/XK**2 

WRI TE ( 6 1 100 1 ) DRCAP, GXO*  DERROR , AO A, DPSERR 

1001  FORMAT ( *  APERTURE  RADIUS3  *, 1PE10.3, 5X, • KAPPA  X3  • ,E10.3,5X, 'ERROR 
1=  • ,E10.3/33X, 'MEAN  SQUARE  ANGLE  OF  ARRIVAL3  • ,E10.3,5X, 'ERROR3  ', 
2E10.3) 

GO  TO  100 
END 
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DOUBLE  PRECISION  FUNCTION  DERFCT (DR ) 

C  LINEAR  PHASE  ERROR  FUNCTION 
IMPLICIT  REAL+8  (D) 

EXTERNAL  DFCT1 
COMMON  R  *RHO, GAMMA* RCA P 
COMMON  / APRAD/DDCAP 
RHC  =DR 

DARG- DR/DDCAP 

CALL  1.MTGRL(1.E-02,1.0E-01*6,DFCT1,DINT1 ) 

DFL=6.*DARCOS(DARG)-(14.*DARG-8.*DARG**3)*  DSQRTI I.-DARG**2 ) 

DERFCT -OR*DFL*OINTl 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  DFCTI ( DXKAP ) 

C  INTEGRAND  FOR  CALCULATING  PHASE  STRUCTURE  FUNCTION  FOR  DERFCT 
IMPLICIT  REAL*8  (D) 

COMMON  R,RHO* GAMMA ,RCAP 

COMMON  /PHI/XKM*XKO 

COMMON  /PHASE/XKtXL 

ARG3* ( OXKAP/XKM) **2 

IF ( ARG3.GT • 174* I GOTO  100 

DARG1»DXKAP**2*XL/XK 

ARG2=DXKAP«'RH0 

IF ( DARGI.EQ.O. )G0  TO  100 

XKAP=DXKAP 

DFCTl= ( l .-BESO t ARG2 ) I  * ( 1 .  +DS IN ( 0ARG1 ) /DARG1 )*PHIN( XKAP ,XKM*XKO) * 
1DXKAP 
RETURN 

100  DFCT1=0. 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  DINTKDR) 
C  R  INTERGRAND 

IMPLICIT  REALMS  (D) 

COMMON  R, RHO*  GAMMA *RCAP 

EXTERNAL  DINT2 

R=DR 

DRCAP*RCAP 

CALL  DQG8  ( 0. DO, DRCAP, DINT2 , DR INT2 ) 
0INTl  =  0RINT2*R*‘R*1W(R! 

RETURN 

END 
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DOUBLE  PRECISION  FUNCTION  DINT2CDRHOI 
C  RHO  INTEGRAND 

IMPLICIT  REALMS  (0) 

COMMON  R»RHOtGAMMA«RCAP 
EXTERNAL  DINT3 
RHO-DRHO 

CALL  DQG4  (O.DO* 3. 1415927D0,DINT3,DRINT3 ) 
0 INT2*DRI NT3*RH0*RH0*W ( RHO ) 

RETURN 

END 


OOUBLE  PRECISION  FUNCTION  D INT 3 ( OGAMMA ) 

C  GAMMA  INTEGRAND 

IMPLICIT  REALMS  (0) 

COMMON  R«RHOtGAMMA«RCAP 
EXTERNAL  DFCT 
GAMMA*DGAMMA 

CALL  I NTGRL I 1 . 0E-03 , l . E-02 • 7# DFC T , DPHASE » 

0 1 NT3*DPHASE*C0S ( GAMMA  I 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  DFCT (DKAP ) 

C  PHASE  STRUCTURE  FUNCTION  INTEGRAND 
IMPLICIT  REAL *6  (0) 

COMMON  R, RHO, GAMMA, RCAP 
COMMON  /PHASE/XK, XL 
COMMON  /PHI /XKM, XKO 

ARG=DKAP*SQRT { R*R+RH0*RH0-2 . *RHO*R*COS ( GAMMA » ) 

XKAP*DKAP 

DARGl»XKAP«2*XL/XK 

DFCT=( I* -BE SOI ARG) )*  PHIN(XKAP,XKM,XKO)*XKAP*l l.+DSIN ( DARGl ) /DARGI 
1) 

RETURN 

END 
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FUNC  f I  UN  PHlN(XKAPfXKMvXKO) 

C  REFRACT I VE  1NOEX  FLUCTAT ION  SPECTRUM 
TEST* ( XKAP/XKM )**2 
IFITEST.GT.174.) GO  TO  100 
A=EXP(-TEST ) 

B*( XKAP**2  +  XK0**2 ) ** 1 11./6.  ) 

TEST l*ALQG10I A I-ALOGIO ( 8 ) 
IF(ABS(TEST1).GT.30.  IGCl  TO  100 
PHIN*.033*A/B 
RETURN 
100  PHIN*0. 

RETURN 

END 


FUNCTION  W(R) 

C  APERTURE  WEIGHTING  FUNCTION 
W*l. 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  DDEN(DR) 
C  DENOMINATOR  INTEGRAND 

IMPLICIT  REAL*8  10) 

R*0R 

DDEN=DR**3*W<R) 

RETURN 

ENO 


SUBROUTINE  lNl  GKUXL$,XUStNINT  ,FCT,YINT) 
C  URUFK  OF  MAGNITUDE  INTEGRATION  ROUTINE 
REAL*b  YINT.YINTG, FCT, XL, XU 
EXTFRNAL  FCT 
XL=XLS 
XU=XUS 

CALL  UCG3? (XLtXUfFCT  »  Y I  NT ) 

UU  300  J=l,NINT 
XL  =  XU 
XU=XU*10. 

CALL  D0G3#>(XLfXUfFCT»YINTG) 
YINT=YINT+YINIG 
300  CONTINUE 
RETURN 
L-NO 


V* 


no 


FUNCTION  BESOUI 
C  JO(X) 

I F ( X.GT. 3* ) GO  TO  100 
Y*X/3. 

B£SP*l.0-2.24999*Y*42+l.26562*Y**4-.3l6387*Y**6+.044447*Y**8-3.944 
14E-34Y4410+2. 10000E-4*Y**12 
RETURN 
100  V*3./X 

THETA=X-. 785398-. 04l6634Y-3.954E-05*Y**2*2.62573E-3*Y**3-5.4l25E-4 
1*Y**4-2.9333E-44Y**5^1.3558E-4*Y4*6 

F 0* .7978846-7 • 7E-7*Y-5. 5274E-3*Y**2- 9.51 2E-5*Y**3*I. 37237E-3*Y**4- 
17.2805E-4*Y665*i.4476E-4*Y**6 
BESO*FO*COS( THETA ) /SORT (X ) 

RETURN 

ENO 
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LISTING  B-6 


C  ****  10.6  MICRON  PROPAGATION  PHENOMENA  PROGRAM  **** 

C  ****  ANCLE  OF  ARRIVAL  ERROR  CALCULATION 

C  USING  THE  PHASE  STRUCTURE  FUNCTION  **♦* 

IMPLICIT  REAL*8  (0) 

COMMON  R, RHO, GAMMA «RCAP * XKAP 
COMMON  /PHI/XKM, XKQ 
COMMON  /PHASE/XK, XL 
EXTERNAL  OFCT 

NAMELIST/ 1 NPUT /XCAPLO, XLQ,RCAP, A, XLAMB«XL 
P I *3. 1415927 
WRITE (6,10001 

1000  FORMAT! 'IANGLE  OF  ARRIVAL  CALCULATION  USING  PHASE  STRUCTURE  FUNCTI 
ION*/*  (RESULTS  NORMALIZED  TO  CN^*2 ) • ) 

100f»  FORMAT!  ‘ORANGE®  *,1PE11.3*'  M',5X, 'INNER  SCALE*  ',E10.3,'  M',5X, 

1 ' OUTER  SCALE*  ',E10.3,'  M',5X,'A=  • , E 10. 3 , 5X, 'LAMBDA*  '.E10.3,'  M« 
2) 

100  RE AD ( 5, I NPUT ) 

WRITF(6,I 005 )XL«XLO«  XCAPLO, A, XL AMB 

XK*6. 28 318/XL AMB 

XKM*» .92/XL0 

XKO-A/ XCAPLO 

DRC AP*KCAP 

CALL  INTGRLI1.0E-03,  l.E-02,7,  DECT ,  DELAOA ) 

DELTAA*  4. *PI **2*XL*DELA0A 
WRITE (6, 1001 »DRCAP,DELTAA 

1001  FORMAT!'  APERTURE  RADIUS*  ' , I PE10. 3, 5X, ' DELTA  AOA*  ',E10.3,/) 

GO  TO  100 

END 


DOUBLE  PRECISION  FUNCTION  DFCT(DKAP) 

C  PHASE  SIRUCTURE  FUNCTION  INTEGRAND 
IMPLICIT  REALMS  (D) 

REAL*8  AUASA* AOAL Av BES IX , AOAL 
COMMON  R « RHO, GAMMA, RCAP, XKAP 
COMMON  /PHASE/XK, XL 
COMMUN  /PHI /XKM, XKO 
EXTERNAL  DINT l 
XKAP=DKAP 

DARG1=XKAP*42*XL/XK 

DRCAP*RCAP 

LALLD0G12 ( 0. DO , DRCAP , D I  NT  1 , DR  I NT l ) 

A0ALA=-16. 40R I  NT  1/ ( 3 • 14 159 2* DRC AP* *8 ) 

AOAL* AOAL A/DKAPP*2 

A0ASA*BESlX(DKAP*DRCAP)/2. 000 00000000000 DO 
DIF*AOASA-AOAL 

UFCT*DIF*UKAP**3*PHIN(XKAP,XKM,XK0)*( I.+DSIN! DARG1 I/DARGl ) 
RETURN 


DOUBLE  PRECISION  FUNCTION  DINT1IDRI 
C  K  INTERGKAND 

IMPLICIT  REAL*8  (0) 

COMMON  R( RHO( GAMMA tRCAPyXKAP 

EXTERNAL  DINI2 

R«DR 

DRCAP*KCAP 

CALL  D0G8  ( O.UO, DRCAP, DINT  2, DRINT2 ) 

DINT1=DRINT2*R*R 

RETURN 

ENO 


DOUBLE  PRECISION  FUNCTION  DINT2I ORHO) 

C  KHQ  INTEGRAND 

IMPLICIT  REAL*8  (0) 

COMMON  R«RHO«GAMMAf RCAPtXKAP 

EXTERNAL  OINT'i 

RHO=DRHO 

CALL  DQG4  ( 0.00, 3. 1M5927D0 , D INT3.UR 1NT3 ) 

U INT2=DRI NT3*RH0*RH0 

RETURN 

END 


DOUBLE  PRECISION  FUNCTION  D INT 3 1 DGAMMA ) 

C  GAMMA  INTEGRANO 

IMPLICIT  REAL*8  (D) 

REAL*8  BESO 

COMMON  R(KHOvGAMMAf RCAPfXKAP 
EXTERNAL  DECT 
GAMMA*OGAMMA 

DPHASE=1 .OOOOOOOOOOOOOOO-BESO ( XKAP^DSQRT ( R**2+RHO**2-2.*R*RHO* 
1DC0S( DGAMMA ) ) ) 

D INT3«DPHASE*C0S (GAMMA ) 

RETURN 

END 


FUNCTION  PHINI  XKAPyXKHt  XKO) 

C  KfcFRACl IVE  INDEX  FLUCTATION  SPECTRUM 
TEST*!  XKAP/XKM ) **2 
IFITESf.GT. 1/4. JGO  TO  100 
A=EXPl -TEST ) 

b»(XKAP**?+XKU**2>**( IL./6. ) 
T£ST1*ALUGI0( A ) - ALUG10 ( b ) 
IFIABSITESTII.GT.30. )G0  TO  100 
PHI N* • 033*A/B 
RETURN 
100  PHI N*0. 

RETURN 

END 


SUBROUTINE  INIGRL (XLStXUStNINTtFGT*Y I NT) 
C  URUEK  OF  MAGNITUDE  INTEGRATION  ROUTINE 
KEAL*fi  Y1NT,Y INTG, FCT ,XL ,XU 
EXTERNAL  FCT 
XL=XLS 
XU=XUS 

CALL  DCG16(XL.XU,FCT,Y1NT) 

DO  300  J= 1 ♦ N I  NT 
XL  =  XU 
XU=XU*10. 

CALL  UCGlto(XL,XU,FCT,YINTG) 
YINT*YINT+YINTG 
300  CONTINUE 
RETURN 
t  NU 


[ 


0 

LI 

D 


L 


FUNCTION  BESO I  X ) 

REAL*B  BESOtX,Y,FO, THETA 
C  JO(X) 

IFIX.GT.i. IGO  TO  100 
Y*X/3. 

BESO- 1  .0-2. 24999* Y**2+l.  26562 *Y* *4-.  3 16387*Y**6+. 04444  7*Y**8-3.94*» 
14E-3*Y**10^2. IOOOOE-4*Y**L2 

return 

100  Y*3./X 

THETA*X-.78539B-.041663*Y-3.954E-05*Y**2*2.625T3E-3*V**3-5.4125E-4 

l#Y**4-2.9333E-4*Y**5+l.3558E"4*Y**6  I 

F0«. 7978846-7. 7E-7*Y-5. 52 74E-3*Y**2-9.5l2E-5*V**34l.37237E-3*Y**4 

17.2806E-44Y**5^l.4476E-4*Y**6 
BESO*FO*OCOS( THETA l/DSQRTI X I 
RETURN 
END 


t 

u 
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FUNGI  ION  BESIX(X) 

REAL*B  THETA. FI. Y.BES1X.X 
C  Ji(X)/X 

IF(X.GT.3. )GO  TU  100 
Y*X/3. 

BESlX».50-.562499*Y**24.?l0935*Y**4-3.95428E-2*Y*4fc*4.433l9E-34Y*4 
18-3.1761E-4*Y**10+1. 109E-5*Y*412 
RETURN 
100  Y-3./X 

THETA«X-2.356l9+.124996*Y+5.650E-05*Y**2-6.3/B?9E-3*Y443+7.434BE-4 

1*Y4*4*7.9824E-4*Y***-2.9166E-4*Y**6 

Fl«. 797884*1. 56E-6*Y  +  i.6!>966E-2*Y**2-l.  71 05E-4*Y4*3-2. 4931 1E-3*Y4* 
14+1.13653E-34Y*45-2.0033E-4*Y**6 
BES IX*F l *DCOS ( THETA )/ ( X*DSQRT IXI) 

RETURN 

END 
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